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Abstract 

Considerable  work  has  been  accomplished  at  AFIT  in  the 
last  three  years  to  improve  the  tracking  capability  of  the 
high  energy  laser  weapon.  The  improvements  were  achieved  via 
use  of  an  adaptive  extended  Kalman  filter  algorithm.  In  this 
research,  work  is  initiated  on  a  tracker  able  to  handle  "mul¬ 
tiple  hot  spot"  targets,  in  which  digital  signal  processing 
is  employed  on  the  FLIR  data  to  identify  the  underlying  target 
shape.  This  identified  shape  is  then  used  in  the  measurement 
model  portion  of  the  filter  as  it  estimates  target  offset  from 
the  center  of  the  f ield-of-view.  Two  tracking  algorithms  are 
developed.  The  first  algorithm  uses  an  extended  Kalman  filter 
to  process  the  intensity  measurements  from  a  FLIR  to  produce 
target  position  estimates.  The  second  algorithm  uses  a  linear 
Kalman  filter  to  process  the  position  estimates  of  an  improved 
correlation  algorithm.  This  algorithm  is  improved  over  stan¬ 
dard  correlators  by  using  thresholding  to  eliminate  poor  cor¬ 
relation  information,  dynamic  information  from  the  Kalman 
filter  and  it  also  uses  the  on-line  derived  target  shape. 
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ENHANCED  TRACKING  OF  AIRBORNE  TARGETS  USING 
FORWARD  LOOKING  INFRA-RED  MEASUREMENTS 

I.  Introduction 

The  first  successful  laser  is  generally  creditted  to 
Maiman  at  the  Hughes  Research  Laboratories  in  1960.  Since 
that  time,  there  has  been  an  explosive  growth  in  laser  tech¬ 
nology.  High  energy  lasers  have  become  prime  candidates  for 
use  as  weapon  systems  (1:87),  and  in  fact  have  already  been 
tested  against  airborne  targets  (2:14-17;  3:16-19) .  The  use 
of  high  energy  lasers  for  military  applications  may  have  a 
significant  impact  on  any  future  war  (1:87). 

Laser  weapon  systems  are  highly  desirable  because  of 
their  potential  to  deposit  large  amounts  of  energy  on  a  small 
area  of  a  target  in  a  short  period  of  time.  However,  there 
are  many  technological  problems  which  must  be  resolved  prior 
to  high  energy  lasers  being  used  as  efficient  and  effective 
weapon  systems. 

The  problem  of  very  accurate  target  position  estimation, 
for  tracking  the  target,  is  one  such  technical  problem.  An 
Infra-Red  (IR)  laser  beam  must  be  kept  positioned  on  a  speci¬ 
fic  part  of  a  target  for  an  extended  amount  of  time  to  destroy 
a  target  or  critical  components  of  a  target.  An  aiming  angle 
correct  to  within  "six  one  hundred  thousandths  of  one  degree", 
.1  microrad,  "is  required  to  hit  a  missile  5000  Km  away,  which 


is  a  precision  beyond  existing  technology  for  both  the  U.S. 
and  the  U.S.S.R."  (1:87-89). 

Background 

This  research  investigates  potential  solutions  to  the 
target  tracking  problems  associated  with  directed  energy  beam 
weapons.  The  Air  Force  Weapons  Laboratory  (AFWL)  located  at 
Kirtland  Air  Force  Base,  New  Mexico,  currently  uses  correlation 
trackers  to  provide  precise  target  position  estimates  to  feed¬ 
back  controllers  in  the  presence  of  several  disturbances. 

These  disturbances  include  any  effect  that  can  cause  relative 
motion  between  the  beam  and  the  target,  such  as  target  motion, 
mirror  vibration,  atmospheric  jitter,  and  sensor  measurement 
errors  (4:2).  A  correlation  tracker  stores  a  complete  set  of 
predetermined  or  previous  real-time  data  which  provide  a  tem¬ 
plate  to  compare  with  new  information  to  be  received  at  a 
later  time.  Cross  correlation  techniques  are  used  to  estimate 
the  relative  position  offsets  from  one  frame  of  data  to  the 
next.  This  relative  position  information  is  then  used  to 
control  the  gimbals  which  keep  the  target  centered  within  the 
sensor's  f ield-of-view.  One  possible  sensor,  the  one  cur¬ 
rently  of  most  interest,  is  the  Forward  Looking  Infra-Red 
sensor  (FLIR) . 

The  high  energy  laser  is  servoed  to  the  FLIR  so  that 
centering  the  target  within  the  FLIR's  f ield-of-view  physi¬ 
cally  points  the  laser  toward  the  target.  This  type  of 
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tracker  needs  no  prior  information  about  the  target’s  shape 
or  dynamics  and  is  therefore  well  suited  to  many  general 
applications.  The  disadvantages  of  this  tracker  are  its  sus¬ 
ceptibility  to  noise  and  its  insensitivity  to  any  knowledge 
of  the  target's  dynamics. 

To  overcome  these  specific  deficiencies  of  the  corre¬ 
lator  tracker,  the  use  of  an  extended  Kalman  filter  tracker 
has  been  investigated  (4; 5; 6).  Appendix  J  contains  a  very 
brief  generic  view  of  the  Kalman  filter  equations  and  then 
provides  information  on  how  they  are  used  for  this  application. 
In  many  practical  tracking  problems,  the  type  of  target  being 
tracked  is  known  along  with  target  parameters  such  as  size, 
shape  and  even  acceleration  characteristics.  Additional  in¬ 
formation  incorporated  by  the  Kalman  filter,  which  is  not 
utilized  by  correlation  trackers,  is  knowledge  of  the  statis¬ 
tical  effects  of  atmospheric  distortion  on  the  radiated  wave- 
front  as  it  propagates  to  the  FLIR.  The  Kalman  filter  uses 
this  information  to  aid  in  separating  the  true  target  relative 
motion  from  the  atmospheric  disturbance.  The  Kalman  filter 
operating  upon  the  raw  digitized  image  has  been  shown  to  per¬ 
form  well  in  comparison  to  correlator  trackers  when  the  tar¬ 
get  intensity  function  shape  is  relatively  well  known  and 
unimodal ,  and  when  the  internal  filter  models  depict  the 
actual  tracking  situation  reasonably  well  (5; 6).  These 
studies  found  the  extended  Kalman  filter  tracker  to  be  a 
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robust  tracker,  even  in  low  signal-to-noise-ratio  environ¬ 
ments  . 

Problem 

A  major  difficulty  in  using  an  extended  Kalman  filter 
for  enhanced  IR  tracking  is  the  Kalman  filter  algorithm's 
requirement  for  an  accurate  reference  image,  h[x(t^),t^]  and 
the  derivatives,  3 [h (x (t^) , t^) ]/3x  =  H(x(t^),t^)  of  that  ref¬ 
erence  image  with  respect  to  the  states,  over  the  entire  FLIR 
image.  This  research  adopts  the  notation  that  an  underlined 
lower  case  letter  denotes  a  vector  while  an  uppercase  letter 
underlined  denotes  a  matrix.  In  previous  research  efforts, 
this  reference  image  was  assumed  to  be  a  known  Gaussian  func¬ 
tional  form.  During  the  past  year,  work  was  initiated  on  a 
tracker  able  to  handle  "multiple  hot  spot"  targets,  in  which 
digital  processing  must  be  employed  to  identify  the  target's 
shape  (4s).  This  identified  shape  could  then  be  used  in  the 
measurement  model  portion  of  the  Kalman  filter  as  it  estimates 
target  offsets  from  the  center  of  the  f ield-of-view.  The 
state  estimates,  in  turn,  are  incorporated  into  the  shape 
function  determination.  The  data  processing  algorithm  is 
shown  in  Figure  1.  The  Kalman  filter  processes  the  measure¬ 
ment  vector  js(t^)  using  the  equation 

x(ti+)  =  x(ti~)  +  K(ti)  [z.(ti)  -hfxU^) ,  ti)  ]  (1-1) 

where 
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x(t .  )  =  state  estimate  vector  after  measure¬ 
ment  incorporation  at  time  t^ 

A  _ 

x(t.  )  =  state  estimate  vector  propagated  from 
previous  measurement  update  to  time  t^ 

K(t^)  =  Kalman  filter  gain 

z.(t.)  =  Measurement  vector  of  average  intensi- 
1  ties  over  individual  picture  elements 
(pixels)  of  the  FLIR  array 

h (x (t . “) , t . )  =  intensity  shape  function  for  measure- 
x  ments  at  the  time  t.  as  a  function  of 
the  state  estimate 

The  apparent  location  of  the  target  within  the  sensor's 
f ield-of-view  is  the  sum  of  effects  due  to  true  target  dyna¬ 
mics  and  atmospheric  disturbances.  The  four-state  estimate 
vector,  xCt^),  consists  of  estimates  for  the  x  and  y  posi¬ 
tions  due  to  true  target  dynamics  and  the  x  and  y  positions 
due  to  atmospherics. 

Figure  1  shows  two  data  processing  paths  for  the  inten¬ 
sity  measurements.  The  FLIR  is  positioned  so  that  the  center 
of  its  field-of-view  is  where  the  Kalman  filter  predicted  the 
true  target  position  to  be,  resulting  from  dynamics  over  the 
most  recent  sample  period.  Each  FLIR  frame  is  arranged  by 
rows  into  a  measurement  vector  z.(t^)  which  is  the  input  to 
the  Kalman  filter  in  the  lower  path  of  Figure  1.  The  Kalman 
filter  uses  the  nonlinear  intensity  function  evaluated  at  the 
predicted  state,  h [x (t^~) , t^J ,  and  the  linearized  intensity 
function,  H[x (t^~) , t^ ] ,  for  that  frame  to  compute  an  updated 
state  estimate,  x  (t^+) ,  from  the  measurement  vector  z.(t^)  . 
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Figure  1.  Data  Processing  Alg< 


The  Kalman  filter  then  uses  its  internal  dynamics  model  to 
propagate  its  state  estimates  forward  to  produce  its  best 
estimate  of  the  states  at  the  next  sample  time,  x(t^~^) • 

That  information  is  passed  to  the  controller,  which  will 
position  the  FLIR  so  that  the  center  of  its  f ield-of-view  is 
where  the  filter  predicts  the  dynamics  will  place  the  target 
position  at  that  next  sample  time.  The  upper  path  of  Figure 
1  is  designed  to  provide  the  linearized  and  nonlinear  inten¬ 
sity  function  for  use  by  the  Kalman  filter. 

The  nonlinear  intensity  pattern,  h  (xft.^-) ,  t^  ,  is  the 
noise-free  intensity  measurement  frame  expected  if  the  Kalman 
filter's  propagated  state  estimates  were  perfect.  Since  the 
center  of  the  f ield-of-view  is  equal  to  the  filter-predicted 
target  location  due  to  the  controller  action,  the  nonlinear 

A  _ 

h[x(t^  ),t^]  will  be  the  target's  intensity  profile,  tj,  off¬ 
set  from  the  center  of  the  FLIR's  f ield-of-view  by  the  pre¬ 
dicted  atmospheric  states.  This  is  because  the  x-axis  posi¬ 
tion  of  the  centroid  of  the  target's  intensity  profile  is 
defined  by 

xcentroid  xdynamics  +  Xatmospherics  (1-2) 

and  the  control  action  is  specifically  intended  to  zero  out 
the  estimated  X£jynamj_cs'  an<^  similarly  for  the  y-axis.  The 

A 

linearized  intensity  pattern,  H [x (ti~) , ti ] ,  is  also  used  by 
the  Kalman  filter,  within  the  computation  of  the  gain  matrix 
K(t^) ,  to  incorporate  a  measurement  to  produce  an  updated 
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state  estimate,  x(t i  ) .  This  linearized  intensity  pattern 
is  the  derivative  of  the  nonlinear  intensity  function, 

A  ^ 

*  t^)  /  with  respect  to  a  change  in  the  Kalman  filter 
states  evaluated  at  the  predicted  state  at  that  sample  time, 
x  (t^~) :  H[x(ti“),ti]  =  3h[x(ti“) , ]/dx. 

To  generate  these  estimated  target  intensity  patterns 
from  the  noise-corrupted  FLIR  requires  interframe  filtering 
to  attenuate  the  noise.  This  inter frame  smoothing  requires 
the  target's  intensity  profile  to  be  centered  from  frame  to 
frame  so  that  the  noise  can  be  averaged  out.  To  center  the 
target's  intensity  function  within  a  data  array  which  repre¬ 
sents  a  FLIR's  field-of-view,  the  shifting  theorem  of  the 
Fourier  Transform  is  used. 

The  offset  of  the  target's  intensity  pattern  from  the 
center  of  the  FLIR's  field-of-view  can  be  negated.  This  is 
accomplished  by  multiplying  the  Fourier  Transform  of  the 
frame  of  data  by  the  complex  conjugate  of  the  linear  phase 
shift  induced  in  the  frequency  domain  by  the  translation  of 
the  intensity  pattern  from  the  center  of  the  frame  in  the 
space  domain.  Atmospheric  jitter  and  imperfect  propagation 
of  dynamic  states  are  the  sources  of  the  image's  translation 
from  the  center  of  the  sensor’s  field-of-view  in  the  space 
domain.  The  output  of  the  Kalman  filter's  update  Equation 
(1-1)  provides  an  estimate  of  both  of  these  offsets  for  the 
present  frame  of  data,  xf^  ) .  These  estimates  of  the  updated 
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target  location  and  atmospheric  jitter,  x (t^+) ,  are  used  to 
compute  the  estimated  centroid  location,  as  shown  by  Equation 
(1-2) .  The  centroid  location  is  used  in  the  argument  of  the 
complex  conjugate  of  the  linear  phase  shift  which  negates 
the  spatial  translation  effects  and  provides  a  centered  in¬ 
tensity  function.  This  centered  pattern  is  then  averaged 
with  previously  centered  frames  of  data,  using  exponential 
smoothing,  to  reduce  the  effects  of  noise  corruption.  The 
derivative  property  of  the  Fourier  Transform  is  then  used  to 
differentiate  the  centered  smoothed  intensity  function  with 
respect  to  a  change  in  Kalman  filter  states. 

This  shape  function  and  its  derivatives  are  then  eval¬ 
uated  at  the  state  expected  at  the  next  sample  time.  For 
this  application,  where  the  FLIR  is  centered  on  the  predicted 
target  location,  the  intensity  patterns  are  evaluated  at  the 
predicted  atmospheric  states.  This  is  the  intensity  pattern 
which  would  be  received  if  the  image  were  noise-free  and  the 
Kalman  filter  had  made  perfect  estimates  of  the  dynamic  and 
atmospheric  states,  of  Equation  (1-2),  which  determine  the 

location  of  the  centroid  of  the  intensity  pattern.  The  in- 

\ 

verse  Fourier  transform  is  then  computed  and  h[x(t^+^) , t^+^] 

*  _  \ 
and  H[x(ti+1) , t^+1 J  are  ready  for  the  Kalman  filter  to  pro¬ 
cess  the  next  frame  of  data. 

Plan  of  Attack 

This  section  provides  an  overview  of  the  general  flow 
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of  this  research,  as  pictured  schematically  in  Figure  2. 


Figure  2.  Research  Plan  of  Attack  Flow  Diagram 


The  first  block  of  Figure  2  represents  the  first  step 
of  this  research.  This  step  includes  the  development  of  the 
algorithms  necessary  to  process  the  FLIR  measurements.  The 
complete  implementation  of  Figure  1  was  accomplished  during 
this  portion  of  the  research.  This  included  the  development 
of  the  four-state  extended  Kalman  filter  needed  to  estimate 
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target  position  along  with  the  estimated  atmospheric  distur¬ 
bance  in  the  x  and  y  directions.  To  generate  the  intensity 
functions  h[x(t^~) , t^]  and  H [x (t^~) , t^ ] ,  algorithms  were 
needed  to  take  the  forward/ inverse  two-dimensional  Discrete 
Fourier  Transform  of  the  FLIR  data.  Once  in  the  frequency 
domain,  the  capability  to  shift  the  data  to  center  it  within 
the  frame  was  developed  so  interframe  filtering  could  be 
accomplished.  An  exponential  smoothing  algorithm  was  devel¬ 
oped  for  filtering  out  background  noise.  The  derivative  of 
this  smoothed  intensity  function  was  derived  via  the  deriva¬ 
tive  property  of  the  Fourier  transform.  The  shifting  algo¬ 
rithm  is  then  used  to  compute  what  the  intensity  function 
and  its  derivative  would  be  as  evaluated  at  the  filter- 
predicted  state  at  the  next  sample  time. 

The  next  objective  of  this  research  was  to  perform  an 
analysis  of  the  target  tracking  (and  shape  function  genera¬ 
tion  secondarily)  performance  and  stability  of  the  data  pro¬ 
cessing  algorithm  of  Figure  1.  Monte  Carlo  analysis  tech¬ 
niques  were  used  to  gather  statistics  on  the  ability  of  this 
data  processing  algorithm  to  generate  h[x(t^-),t^]  and 
H [x. ( ) , t^] .  It  is  the  target  tracking  performance  which  is 
of  primary  concern. 

The  alternate  data  processing  algorithm  evaluation 
section  evaluates  the  potential  of  an  improved  correlator 
tracker  followed  by  a  Kalman  filter.  The  correlator  tracker 
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works  well  in  high  signal-to-noise  ratio  scenarios  and  pro¬ 
vides  good  estimates  of  centroid  location,  especially  for 
extended  targets  (as  opposed  to  point  source  targets) .  The 
Kalman  filter  could  then  filter  these  estimates,  treated  as 
measurements  for  the  filter,  to  separate  out  the  components 
of  Equation  1-2. 

The  demand  for  real-time  operation  of  the  computationally 
intense  algorithms  discussed  above  requires  the  development 
of  parallel  processing  techniques  such  as  optical  processing. 
The  final  section  is  a  study  into  the  potential  of  optically 
computing  some  of  the  quantities  of  Figure  1. 

Overview 

Chapters  II,  III,  and  IV  will  describe  in  detail  the 
mathematical  models  used  in  the  computer  implementation. 

Chapter  II  presents  the  algorithms  used  in  the  derivation  of 
the  nonlinear  and  linearized  intensity  functions.  Chapter 
III  presents  the  mathematical  model  used  for  the  truth  model 
which  represents  the  environment  from  which  measurements  are 
taken.  Chapter  IV  describes  the  Kalman  filter  model  used  in 
the  computer  simulation.  Chapter  V  discusses  the  possibility 
of  a  correlator  followed  by  a  Kalman  filter  to  provide  the 
position  estimates.  Chapter  VI  presents  the  results  from 
the  performance  analysis  which  includes  the  statistics  on 
how  well  the  intensity  functions  are  being  constructed  as 
well  as  statistics  on  target  tracking  ability  for  both 
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algorithms#  and  Chapter  VII  presents  the  optical  processing 
options.  Finally,  Chapter  VIII  will  present  the  conclusions 
and  recommendations. 
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II.  Derivation  of  the  Nonlinear  and 
Linearized  Intensity  Functions 

Introduction 

This  chapter  will  present  the  algorithms  necessary  to 
produce  an  accurate  reference  intensity  profile,  h [x (ti_) , ^ ], 
and  the  derivatives  of  that  function,  H[x (ti~) , t^ ] ,  with  re¬ 
spect  to  the  states.  These  intensity  functions  are  needed  by 
the  extended  Kalman  filter  to  process  FLIR  measurements  to 
update  state  estimates  (see  Chapter  IV,  Extended  Kalman  Fil¬ 
ter)  .  The  true  intensity  profile  is  a  vector  of  scalar  com¬ 
ponents  equal  to  the  average  intensity  of  the  target's  IR 
image  over  a  particular  pixel  in  the  f ield-of-view.  Noise- 
corrupted  FLIR  data  must  be  processed  to  determine  the  tar¬ 
get's  profile.  Figure  1  showed  the  data  processing  necessary 
to  produce  the  intensity  functions  from  the  incoming  data. 

The  goal  is  to  recognize  the  intensity  function  which  is 
common  to  the  noise-corrupted  frames  of  data. 


Pattern  Recognition 

All  of  the  information  of  a  two-dimensional  intensity 

/ 

pattern  can  be  preserved  by  a  set  of  eigenvalues  and  eigen¬ 
functions.  To  retain  all  of  the  information  of  a  profile 
may  require  an  infinite  set  of  such  values  and  functions. 

It  is  desirable  to  examine  the  FLIR  image  data  in  a  coordinate 
system  which  has  properties  more  conducive  to  recognizing 
patterns  than  the  spatial  x-y  coordinate  system.  Ideally,  a 


linear  transformation  which  rotates  the  FLIR  image  into  a 

new  coordinate  space  where  the  components  are  uncorrelated 

is  desired.  The  Karhunen-Loeve  Transformation  is  an  example 

of  such  an  operation.  This  transformation  does  possess  the 

disadvantage  of  having  to  produce  the  correlation  matrix 
2  2 

which  is  N  x  N  if  the  input  square  matrix  is  N  x  N. 

The  Karhunen-Loeve  Transformation  possesses  the  pro¬ 
perty  of  providing  a  coordinate  space  with  perfectly  uncor¬ 
related  image  components,  but  it  is  a  difficult  transformation 
to  perform  in  its  exact  form.  The  Karhunen-Loeve  Transforma¬ 
tion  does  provide  motivation  for  the  use  of  the  more  familiar 
Fourier  Transform  which  uses  complex  exponentials  as  eigen¬ 
functions  (4:12).  Where,  in  the  Karhunen-Loeve  Transformation 
the  eigenvalues  correspond  to  actual  variance  statistics  pro¬ 
jected  on  coordinate  axes  of  basis  vectors  in  a  space  where 
the  image  data  is  uncorrelated,  the  Fourier  Transform  is  a 
projection  of  the  image  along  the  basis  vectors  formed  by 
complex  exponentials  (9:195). 

Although  the  Fourier  Transform  does  not  provide  the 
perfect  decorrelation,  it  does  possess  a  computational  advan¬ 
tage  which  outweighs  the  compression  efficiency  of  the 
Karhunen-Loeve  Transformation  (9:194).  The  Fourier  Transform 
also  possesses  the  advantage  that  the  two-dimensional  trans¬ 
form  is  separable  and  can  be  obtained  by  one-dimensional 
operations. 
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Two-Dimensional  Fourier  Transforms 

The  purpose  of  decomposing  the  information  contained 

within  the  complicated  image  intensity  data  into  a  set  of 

eigenvalues  and  eigenfunctions  is  to  compress  that  information 

into  a  set  of  inputs  which  will  make  manipulations  performed 

on  the  data  easier.  Fourier  analysis  provides  one  way  of 

performing  such  a  decomposition. 

The  Fourier  Transform  of  a  complex-valued  function 

g(x,y)  of  two  independent  variables  is  a  decomposition  into 

a  linear  combination  of  elementary  functions  of  the  complex 

exponential  form  exp[j2Tr(f  x  +  f  y)  ]  (7:8).  This  Fourier 

x  y 

Transform  is  defined  by 

G  ( f  ,  f  )  =F  ( g  (x ,  y ) )  =fj  g  (x,y)  exp[-j2TT(f  x  +  f  y)  ]dxdy  (2-1) 
x  y  *  y 

where 

G (f  ,  f  )  =  Frequency  Spectrum 
x  y 

g(x,y)  =  Spatial  Domain  Representation 

f  ,f  =  Spatial  Frequencies 
x  y 

x,y  =  Spatial  Variables 

The  transform  is  also  a  complex-valued  function  of  two 
independent  variables  fx  and  fy/  which  are  the  spatial  fre¬ 
quencies.  The  inverse  Fourier  Transform  is  defined  by 

g(x,y)=F-1(G(f j/y)  )=//  G(fx,fy)exp[+j27i(fxx  +  fyy)  ]dfxdfy 

(2-2) 

The  complex- valued  number  G(f  , f  )  is  a  weighting  factor  that 

x  y 
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is  applied  to  the  eigenfunction,  which  represents  the  coor¬ 
dinate  axes  of  the  Fourier  domain,  in  order  to  synthesize 
g(x,y) . 

The  Fourier  Transform  is  separable,  so  that  the  two- 
dimensional  transform  domain  can  be  reached  with  one-dimen¬ 
sional  operations.  Implementation  is  accomplished  by  trans¬ 
forming  first  on  the  rows  of  the  image,  followed  by  the  trans¬ 
formation  along  the  columns. 

The  two-dimensional  finite  Discrete  Fourier  Transform 
of  a  two-dimensional  periodic  sequence,  of  period  N  in  both 
directions,  is  the  finite  sum  of  complex  exponentials  in  the 
form 

G(fx,f  )=  I  Z  g(x,y)expC-j“jL(xfx  +  yf  )  ]  (2-3) 

*  x=0  y=0  y 

The  complex  sequence  of  intensity  values,  g(x,y)  is  discre¬ 
tized  into  an  N-pixel  by  N-pixel  array.  The  inverse  Discrete 
Fourier  Transform  is  defined  by 

N-l  N— 1  „  0tt 

gU,y)=-4-  I  Z  G(f  ,f  )exp[+#(f  x  +  f  y)]  (2-4) 

Nf  =0  f  =0  X  y  wx  y 

x  y 

The  two-dimensional  complex- valued  N  x  N  array,  G(fx, fy), 
which  represents  the  two-dimensional  finite  Discrete  Fourier 
Transform  of  the  original  FLIR  data,  is  discretized  into 
spatial  frequency  bins.  Each  component  of  the  array  represents 
one  of  the  N/2  spatial  frequencies  or  their  conjugates.  The 
inverse  transform  differs  from  the  forward  transform  only  in 
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the  sign  of  the  exponent  which  appears  in  the  complex  expon¬ 
ential  within  the  summation  of  Equation  (2-4) .  The  transform, 

G(f  ,f  ),  will  be  periodic  of  period  N  as  was  the  original 
x  y 

sequence,  g(x,  y).  Appendix  A  contains  a  further  interpreta¬ 
tion  of  the  two-dimensional  Discrete  Fourier  Transform.  It 
also  contains  details  on  how  the  two-dimensional  Fourier 
Transform  implementation  was  tested. 

In  summary,  the  two-dimensional  Discrete  Fourier  Trans¬ 
form  assumes  a  finite  length  two-dimensional  sequence  is  one 
period  of  a  two-dimensional  periodic  sequence  which  then  can 
be  decomposed  into  a  linear  combination  of  complex  exponen¬ 
tials.  See  Appendix  A  for  explicit  details. 

Shifting  Property  of  Fourier  Transforms 

To  generate  the  target  intensity  patterns  from  the  noise- 
corrupted  FLIR  data,  interframe  filtering  is  required  to 
attenuate  the  noise.  Interframe  smoothing  requires  that  the 
target's  intensity  profile  be  centered  from  frame  to  frame. 
Successive  centered  frames  of  data  can  then  be  averaged  to 
attenuate  the  corrupting  noise.  The  shifting  property  of  the 
Fourier  Transform  is  used  in  conjunction  with  the  filter's 
estimated  location  of  the  intensity  profile  to  center  the 
data. 

A  shift  in  the  scalar  spatial  domain  of  a  finite  dura¬ 
tion  one-dimensional  sequence  can  be  interpreted  as  a  rotation 
in  the  basic  interval  (9:103).  That  is,  the  shifting  of 


18 


mmtwrmxsmiP 


samples  out  of  one  side  of  the  periodic  sequence  results  in 
identical  samples  entering  the  interval  at  the  other  end. 

This  is  sometimes  called  a  cylindrical  shift.  if,  for  example, 
the  Fourier  Transform  of  g(x)  is  defined  by 

F[g(x)  ]=G(u)  (2-5) 

then  the  Fourier  Transform  of  the  shifted  sequence  is  defined 

by 

F[g(x-xo)  ]=e"ja,xo  G(co)  (2-6) 

In  the  case  of  a  finite-area  two-dimensional  sequence,  a 

translation  of  the  function  in  the  space  domain  introduces  a 
linear  phase  shift  in  the  frequency  domain  (9:9): 

F[g(x-xo,y-yQ)  ]=G(fx,fy)exp(-j2ir(fxxo  +  fyyQ)  (2-7) 

where  Fig (x, y) ]=G (f  , f  ) 

a  y 

An  interpretation  of  the  translation  by  xQ  and  yQ  in 
the  space  domain,  analogous  to  the  one-dimensional  interpre¬ 
tation,  would  be  a  rotation  of  each  column  by  xq  samples  and 
each  row  by  yq  samples  (9:119). 

The  offset  of  the  target's  intensity  pattern  from  the 
center  of  the  FLIR's  f ield-of-view  can  be  negated  to  provide 
a  centered  image  for  smoothing.  This  is  accomplished  by 
multiplying  the  Fourier  Transform  of  the  frame  of  data,  which 
contains  the  offset  intensity  profile,  by  the  complex  conju¬ 
gate  of  the  linear  phase  shift  induced  in  the  frequency 
domain.  The  centered  intensity  profile  can  be  obtained  from 
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the  offset  profile  by 


g(x/y)=F-1{F[g(x-x  «y-y  )  Jexp[+j27r(f  x  +  f  y  )]}  (2-8) 

u  u  a  u  y  u 

In  summary ,  by  using  the  Kalman  filter's  updated  esti¬ 
mates  of  the  location  of  the  centroid  within  a  frame  of  data, 
the  centered  intensity  profile  can  be  obtained.  Appendix  B 
contains  the  details  of  how  the  shift  is  implemented. 

Smoothing  in  the  Fourier  Domain 

The  target's  intensity  profile  is  not  directly  available 
for  observation.  This  is  because  each  FLIR  frame  is  corrupted 
by  inherent  FLIR  errors  and  background  noise.  Background 
noise  can  vary  from  zero  to  high  temporal  and  spatial  correl¬ 
ations.  The  FLIR  errors,  such  as  thermal  noise  and  dark  cur¬ 
rent  effects,  can  be  modelled  as  temporally  and  spatially 
uncorrelated  noise.  The  effects  of  these  noises  on  target 
intensity  profile  shape  generation  must  be  minimized  by  some 
appropriate  data  processing  techniques.  Any  chosen  data  pro¬ 
cessing  technique  must  not  assume  the  precise  form  of  the 
corrupting  noise  because  of  the  large  range  of  correlation 
characteristics.  However,  it  is  appropriate  to  exploit  the 
fact  that  the  corrupting  noise  changes  faster  from  sample 
period  to  sample  period  than  the  target's  intensity  pattern 
(4; 5)  . 

An  exponential  smoothing  algorithm  was  chosen  to  combat 
the  effects  of  the  corrupting  noise.  The  equation  which 
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defines  exponential  smoothing  is: 

£(t)=ay(t)  +  (l-a)y(t-l)  (10:115)  (2-9) 

where 

A 

y(t)  =  current  averaged  value 

y (t)  =  current  data  frame 

y(t-l)  =  previous  averaged  data  frame 
a  =  smoothing  constant  o£a£l 
The  smoothing  constant  a  can  vary  anywhere  from  0  to  1. 
The  value  of  a  determines  how  the  current  averaged  data  frame, 
y(t),  responds  to  the  current  data  frame,  y(t).  If  the  tar¬ 
get's  intensity  pattern  is  only  changing  slowly,  the  value  of 
a  should  tend  more  toward  zero  than  one  so  that  more  emphasis 

is  placed  on  the  previous  averaged  data  but  the  present  frame 

is  still  used. 

A 

In  summary/  y (t)  contains  information  from  all  the  past 
data  frames.  The  initial  data  frames  have  less  of  an  effect 

A 

on  the  most  recent  averaged  data  frame,  y(t),  than  do  the  most 
recent  data,  as  is  appropriate  for  a  slowly  changing  target 
image.  Appendix  C  contains  the  details  of  how  this  smoothing 
algorithm  was  implemented  and  how  the  algorithm  can  be  ex¬ 
pressed  in  a  closed  form. 

Derivative  Property  of  Fourier  Transforms 

As  stated  in  Chapter  I  and  to  be  shown  explicitly  in 
Chapter  IV,  the  extended  Kalman  filter  algorithm  requires 
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the  derivatives  of  the  intensity  function  with  respect  to 
the  states,  H[x  (ti_) ,  t^  =  3h  (x  (t^) ,  t^/Sx.  In  past  re¬ 
search  efforts  (4:),  this  linearized  intensity  function  was 
derived  using  a  numerical  approximation  known  as  the  Forward- 
Backward  Difference  Method.  This  method  only  uses  the  pixels 
just  before  and  just  after  the  pixel  being  computed  in  the 
direction  the  spatial  derivative  is  being  taken,  it  was  for 
this  reason  that  the  Derivative  Property  of  the  Fourier  Trans¬ 
form,  which  uses  all  the  data  array,  was  chosen  as  a  better 
means  of  generation  for  this  research. 

Differentiation  in  the  space  domain  just  becomes  a  mul¬ 
tiplication  by  j2ir(f  +  f  )  in  the  frequency  domain.  The 

x  y 

derivative  of  the  intensity  function  can  be  expressed  in 
terms  of  the  transform  of  the  function 

=  j2irfx  •  F[h(x,y)]  (2-10a) 

=  j2irfy  •  F  [h  (x,  y)  ]  (2-10b) 

Appendix  D  provides  explicit  details  on  the  implementation 
of  this  algorithm. 

In  summary,  the  derivatives  of  the  intensity  function 
with  respect  to  the  states  are  required  by  the  extended  Kalman 
filter  algorithm.  These  derivatives  can  be  expressed  in  terms 
of  the  transform  of  that  intensity  function  using  Equations 
(2-10a)  and  (2-10b) .  Appendix  D  contains  the  details  of 


3h (x,y) 
3x 

~3h  (x,  yj 
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Subroutine  Deriv. 


Derivation  of  Intensity  Functions:  Summary 

This  section  will  review  the  algorithms  used  to  produce 
an  accurate  reference  intensity  profile,  h [x (t^-) , t^ ] ,  and 
the  derivatives  of  that  function  H [x (ti~) , t^ ] .  As  stated  in 
Chapter  I,  the  extended  Kalman  filter  algorithm  requires 
these  intensity  functions  for  the  update  state  equations, 

(see  Chapter  IV,  Extended  Kalman  Filter) . 

Noise-corrupted  FLIR  data  must  be  processed  to  deter¬ 
mine  the  target's  intensity  profile.  The  upper  path  of  Figure 
1  showed  the  data  processing  necessary  to  produce  the  inten¬ 
sity  functions  from  the  incoming  noise-corrupted  FLIR  data. 

The  goal  of  the  data  processing  is  to  recognize  the  intensity 
function  which  is  common  to  the  noise-corrupted  frames  of 
data. 

Interframe  smoothing  is  required  to  generate  the  esti¬ 
mated  target  intensity  pattern  from  the  FLIR  data.  This 
smoothing  is  accomplished  in  the  frequency  domain  since  that 
is  where  the  shift  is  accomplished,  so  the  two-dimensional 
transform  is  first  taken  of  the  data.  It  will  be  shown  in 
Appendix  C  on  exponential  smoothing  that  smoothing  in  the 
space  domain  is  equivalent  to  smoothing  in  the  frequency 
domain  (see  Table  I,  Chapter  VI) .  To  average  from  frame  to 
frame,  the  target's  intensity  profile  must  be  centered.  This 
is  accomplished  via  use  of  the  Shifting  Theorem  of  the 
Fourier  Transform.  Once  the  data  is  centered,  exponential 
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smoothing  is  used  to  attenuate  the  noise  while  still  in  the 
frequency  domain.  The  spatial  derivatives  are  then  obtained 
from  this  smoothed  data  by  using  the  derivative  property  of 
the  Fourier  Transform.  The  shifting  property  can  then  be 
used  again  to  derive  the  intensity  functions  evaluated  at  the 
predicted  states.  The  high  energy  laser  will  be  pointed  at 
the  best  estimate  of  the  target's  position  as  calculated  by 
the  Kalman  filter.  The  target's  intensity  function  will  be 
offset  by  atmospheric  jitter  from  that  target  location.  Since 
the  field-of-view  is  centered  at  the  estimated  target  location, 
the  algorithm  expects  the  intensity  function  to  be  offset  from 
the  center  of  the  field-of-view  by  the  estimated  atmospheric 
jitter.  By  using  the  shifting  property  the  estimated  centered 
target  intensity  functions  can  be  manipulated  to  produce  the 
expected  offset  intensity  functions. 

In  summary,  this  chapter  presented  the  algorithms  used 
to  produce  the  intensity  functions  needed  by  the  extended 
Kalman  filter  algorithm.  The  use  of  this  procedure  eliminates 
the  need  to  make  assumptions  on  the  target's  intensity  profile. 
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III.  Truth  Model 


Introduction 

The  truth  model  is  the  best  mathematical  representation 
of  the  real  world  environment  from  which  the  measurement  vec¬ 
tor,  .zft^ ,  is  generated.  The  environmental  characteristics 
which  are  modelled  include  the  underlying  target  dynamics, 
the  atmospheric  jitter  effects,  and  the  background  and  FLIR 
noises.  This  chapter  presents  the  models  used,  along  with  an 
explanation  of  how  the  measurement  vector,  z. (t^) ,  was  created. 
The  first  two  sections  of  this  chapter  present  the  two  pro¬ 
cesses  which  contribute  to  movement  of  the  intensity  function, 
target  dynamics  and  atmospheric  jitter.  These  models  are  then 
transformed  into  state  space  notation  and  the  propagation 
equations  are  given.  Information  on  the  computer  creation 
of  the  measurement  data  for  the  computer  simulation  is  next, 
with  the  final  section  containing  the  model  for  spatial  cor¬ 
relations  of  background  noise. 

Target  Dynamics  Model 

For  this  research,  both  deterministic  and  stochastic 
dynamic  models  were  used.  In  the  deterministic  model,  a 
target  moving  across  the  image  plane  at  a  constant  velocity 
was  simulated.  In  the  stochastic  model,  a  first  order  Gauss- 
Markov  process  was  used  to  describe  the  movement  of  the  tar¬ 
get  with  respect  to  the  f ield-of-view  of  the  sensor.  Although 
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there  exist  better  models  for  describing  the  movement  of 
specific  targets,  the  very  general  stochastic  model  is  appli¬ 
cable  to  many  general  targets.  Together,  these  two  models 
yield  some  of  the  basic  characteristics  of  trajectories. 

The  deterministic  model  was  developed  to  test  the  algo¬ 
rithm  of  Figure  1  with  a  target  which  is  maintaining  a  con¬ 
stant  velocity  across  the  FLIR  sensor's  image  plane.  In  this 
project,  a  constant  velocity  of  three  pixels  over  any  twenty- 
frame  time  history  was  used. 


(ti+l}  =  xd(ti)  +  *15 

(W  =  yd(V  +  *15 


(3-1) 


The  stochastic  model  developed  in  Mercier's  thesis, 
(5:9-10),  for  target  dynamics  was  also  used.  In  this  model, 
the  output  of  a  first  order  linear  shaping  filter  driven  by 
zero  mean  white  Gaussian  noise  was  used  to  describe  the  move¬ 
ment  of  the  target  with  respect  to  the  FLIR's  field-of-view. 
The  outputs  of  these  first  order  shaping  filters,  which  de¬ 
scribe  the  target  dynamics  in  each  direction,  can  be  expressed 
as 

xd(t)  =  ^  *d(t)  +  ^(t)  (3-2) 

yd(t)  =  T^  Yd(t)  +  W2(t)  (3"3) 

where 

E[w1(t)]  =  E[w2(t)]  =  0  (3-4) 
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2 


(3-5) 


Efw^tjw^s)  )  =  E[w2  (t)w2  (s)  ]  = 


6 (t-s) 


EtWj^  (t)w2  (s)  ]  =  0 


(3-6) 


and 


Tfc  =  target  correlation  time 

w, (t) ,w?(t)  =  continuous  time,  independent,  white 
Gaussian  noise  processes 

2 

0^  =  desired  variance  on  outputs  and  y^ 


Atmospheric  Jitter  Model 

This  section  presents  the  mathematical  model  of  the 
translational  position  changes  of  the  target's  intensity  dis¬ 
tribution  due  to  disturbances  in  the  atmosphere  causing  phase 
distortions  in  the  radiated  wavefronts.  The  model  used  is 
based  on  a  study  accomplished  by  The  Analytic  Sciences  Cor¬ 
poration  (11:29-30)  and  data  analysis  by  Hogge  and  Butts 
(12:).  As  a  part  of  these  studies,  a  third  order  shaping 
filter  driven  by  white  Gaussian  noise  was  developed  to  gener¬ 
ate  an  output  with  a  power  spectral  density  representation 
that  well  approximated  that  of  the  effects  of  atmospheric 
turbulence.  The  jitter  of  the  target  intensity  in  each  in¬ 
dependent  direction  is  modelled  as  the  output  of  the  third 
order  shaping  filter  shown  in  Figure  3  (5:12) .  This  filter 
has  a  single  pole  at  14.14  rad/sec  and  a  double  pole  at  659.5 
rad/sec,  and  it  is  driven  by  white  Gaussian  noise  that  is 
zero  mean  and  has  a  strength  of  unity: 
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w  _  > 

2 

K“lw2 

Y 

(s+w1)  (s+0)2)  2 

■  -  * - -  ■? 

where 


w  =  zero  mean,  unit  strength,  white  Gaussian 
noise  process 

,  to_  =  break  frequencies  (u>,  =  14.14  rad/sec, 

1  *  w2  =  659.5  rad/sec) 

y  =  output  of  system 
K  =  system  gain 


Figure  3.  Third  Order  Shaping  Filter 
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E [w  (t) ]  =  0  (3-7) 

E[w(t)  w(s) ]  =  16  (t-s)  (3-8) 

As  shown  by  Mercier,  (5),  by  adjusting  the  value  of  the  system 
gain,  K,  the  desired  root  mean  squared  (RMS)  atmospheric 
jitter  characteristic  is  achieved.  Appendix  C  of  Mercier 's 
thesis  contains  the  details  of  this  derivation. 

State  Space  Model 

The  mathematical  models  developed  in  this  chapter  will 
now  be  transformed  into  state  space  notation.  This  transfor¬ 
mation  results  in  a  time  invariant  vector  differential  equa¬ 
tion  of  the  form 

x^t)  =  F^U)  +  G^U)  (3-9) 

where 

Fl  =  truth  model  plant  matrix 
x^t)  =  truth  model  state  vector 

=  truth  model  noise  input  matrix 
w^  =  vector  of  white  Gaussian  noise  inputs 

E  (Wj.  (t)  ]  =  0  (3-10) 

Etw^tjwjtt  +  t)]  =  ^.6(1)  (3-11) 

The  truth  model  state  vector  is  defined  by 
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=  14.14  rad/sec  =  659.5  rad/sec  (3 


The  truth  model  noise  matrix  (5:11-14)  is  defined  by 


(uij-UjH  (w1-w2) 

Using  Equations  (3-5)  and  (3-8) ,  of  Equation  (3-11) 
becomes 


(3-15) 


The  movement  of  the  target's  intensity  profile,  the  centroid, 
is  governed  by  components  of  these  truth  model  matrices 
where 


xpeak 


(t)  =  x,  (t)  +  x  (t) 

d  a 


=  Y/.(t)  +  y  (t) 


(3-16a) 


xd(t)  =  xt(l) 

xa(t)  =  xt(2)  +  xt(3) 
yd(t)  =  xt(5) 

Ya(fc)  =  xt (6)  +  xfc(7)  (3-16b) 


Thus,  the  location  of  the  centroid  of  the  intensity  profile 
can  be  expressed  as 


yt  (t)  = 

xpeak(t) 

= 

- 

1  1 

1  0  0  0  0  0 

.*peak(t)_ 

0  0 

0  0  0  1  1  1 

(3-17) 

Propagation  Equations 

Propagation  equations  describe  how  the  state  vector 
makes  the  transition  from  one  sample  time  to  another.  The 
solution  to  Equation  (3-9)  (5:14)  is 


fci+l. 


*t<w  =  it'WVW  \!  W i-T)St(T)Kt(T)dT 

(3-18) 


where  _£t(t,t^)  is  the  solution  to  the  matrix  differential 
equation 

-t(t,ti}  =  4it(t,ti!  (3-19a) 

and  initial  condition 

— t^i^i^  =  —  (identity  matrix)  (3-19b) 

The  matrix  $^.(t^+^,t^)  is  t^ie  sta^e  transition  matrix  which 
describes  the  truth  model  state's  transition  from  time  t^  to 
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time  tj_+ji  •  The  truth  model  plant  matrix,  F^.  of  Equation 
(3-14),  is  a  constant  matrix,  and  therefore,  the  state  tran¬ 
sition  matrix  is  only  a  function  of  the  time  difference 
(t..,-t,).  Therefore,  the  state  transition  matrix  becomes 
a  constant  for  a  fixed  sample  time  At  =  ti+1~t^.  The  inde¬ 
pendent  channels  reflected  in  the  truth  model  plant  matrix 
result  in  the  state  transition  matrix  being  block  diagonal  with 
both  of  the  4-by-4  diagonal  blocks  equal  to 


V 


(ti+l' 


V 


e-At/Tt 


0 

0 


0 

-u>,  At 
e  1 

0 


0 

0 


e-u)2At 


0 


Ate 

e 


0 

0 

-0)2At 

— u2At_ 

(3-20) 


where  At  =  constant  sampling  time  (t^+^-t^)  •  A  probabilistic 
interpretation  of  the  integral  term  of  Equation  (3-18)  is 
required  to  describe  the  process  x^.(t).  This  integral  has  a 
zero  mean  Gaussian  distribution.  For  notational  convenience 
the  contribution  of  the  integral  is  set  equal  to  w^(t^)  . 


-i+1 
/ 
t 


it(ti+i'T)^ 


(r) w,  (t) dr 


(3-21) 


Computing  the  statistics  of  v^^Ct^)  (5:15)  results  in 

Etw^u.n  =  o 


(3-22) 
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fci+l 

Elstdlti)^td(ti)1  =  [  £t(ti+i'T)^t(T)Qt(x)^t(T)it(ti+i'T)dT 

r 

(3-23) 
(3-24) 


The  equivalent  discrete-time  model  (5:15)  for  the  above 
equations  is  given  by  Mercier  to  be 


*t<*W  =  +  9gtd  VV 


(3-25) 


where  /2t(}  is  the  Cholesky  square  root  (13:  370)  of 


A  i+1  T  T 

Qtd"  1  £t(ti+l'T)s±(T^t(T^(T)£t(ti+i'T)dT 
fci 


(3-26) 


and 


E[w  (t. ) ]  =  0 
— n  i  — 


(3-27) 


'  I5ij 


(3-28) 


The  lower  triangular  Cholesky  square  root  of  satisfies 


the  equation 


^=d  hi  -  std 


(3-29) 


and  can  be  uniquely  defined  for  a  given  Q 


•td' 


Measurements 

A  measurement  history  of  FLIR  data  frames  is  used  by 
the  Kalman  filter  to  predict  the  pointing  direction  for  the 
high  energy  laser.  The  output  of  an  arbitrary  pixel,  say 
the  k-£-th  pixel,  within  the  FLIR  data  frame  is  the  average 
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IR  intensity  oyer  that  pixel  as  sensed  by  a  detector,  and  is 
described  mathematically  as 


kf/.thItarget(z-5'-xpaak<ti,'lrpeak(tt,,dxdir  +  Vkl<V 

pixel 


(3-30a) 


where 


target 


=  the  intensity  of  the  target  as  a  function  of 
the  location  in  the  frame  and  the  centroid 


peak  as  defined  by  Equation  (1-2) 

(t±)  =  the  output  of  the  k-Jt-th  pixel  at  a  time  t^ 

A  =  the  area  of  one  pixel 
P 

(x,y)  =  coordinates  of  any  point  in  the  k-&-th  pixel 
(t.),y  ,(t.))  =  the  center  of  the  Gaussian  in- 


(x 


peak  i  ' 1  peak  i 


tensity  pattern 


vkj,(t.)  =  the  additive  noise  for  the  k-fc-th  pixel  cor- 
x  responding  to  background  and  FLIR  noises. 


in  the  general  case.  Using  a  bivariate  Gaussian  target  in¬ 
tensity  profile  with  circular  equal  intensity  contours  as  a 
special  case  (4:25),  the  mathematical  description  becomes 


•M"l'  -  Ap  k:e:th  W,xP(!?g  '<x-xpeak(ti))2 


(t  J  =  t  IS 

k-Z-i 
pixel 


+  (y-ypeak(ti))2lldxdy  +  vkil(ti) 


(3-30b) 


where 

I  =  the  maximum  intensity  received  from  target 
max  J 

2 

a g  =  the  dispersion  of  the  Gaussian  intensity  function 
This  would  be  the  measurement  vector  for  a  single  Gaussian 
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hotspot  target.  For  simulation  of  the  multiple  hotspot  tar¬ 
gets  of  interest  for  this  research,  three  exponentials  of 
Gaussian  form  are  summed  in  place  of  the  single  Gaussian  of 
Equation  (3-30a) . 

Spatially  Correlated  Background  Noise 

The  existence  of  spatial  correlations  in  the  background 
of  each  FLIR  data  frame  has  been  documented  by  Harnly  and 
Jensen  (6:14).  To  discuss  the  generation  of  the  spatially 
correlated  noise,  a  numbering  system  for  each  pixel  of  an 
8x8  input  array  is  adopted  (Figure  4).  The  8x8  FLIR 
data  array  is  numbered,  as  done  by  Harnly  and  Jensen  (6:19) 
and  Singletery  (4:21)  by  pixels  from  1  to  64  starting  in  the 
upper  left-hand  corner  and  proceeding  across  the  rows. 

In  the  64  x  64  correlation  coefficient  matrix  corre¬ 
sponding  to  a  noise  vector  with  each  component  associated 
with  a  single  such  pixel,  each  element  relates  the  noise  in 
that  pixel  with  the  noise  in  another  pixel.  For  this  re¬ 
search  the  spatial  correlations  are  modelled  as  extending 
to  the  first  and  second  nearest  pixels.  The  element  in  row 
1  and  column  1  of  the  correlation  coefficient  matrix  relates 
the  noise  in  pixel  one  to  itself  and  is  therefore  one. 
Similarly,  the  entire  diagonal  of  the  correlation  coefficient 
matrix  is  equal  to  one.  The  correlation  coefficient  matrix 
is  symmetric  since  the  correlation  r^j  of  the  noise  in  pixel 
i  with  noise  of  pixel  j  is  the  same  as  that  of  the  noise  in 
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Figure  4.  FLIR  Data  Frame  Pixel  Numbering  Scheme  (6:19) 

pixel  j  with  the  noise  in  pixel  i  (r.^.  =  r^)  . 

The  exponentially  decaying  and  radially  symmetric  form 
of  the  correlation  function  was  used  in  generating  the  coef¬ 
ficients  of  the  first  and  second  nearest  neighbor  correlation 
matrix.  Equation  (3-31) ,  was  used  to  generate  the  coefficient 
matrix: 

r.  .  =  exp(-d.  .)  (3-31) 

ij  *  ij 

where  d^j  =  distance  in  #  pixels  from  pixel  i  to  pixel  j . 
Equation  (3-31)  implements  a  correlation  distance  of  one 
pixel.  The  fact  that  nonzero  values  are  computed  only  for 
the  first  and  second  neighbors  results  in  only  the  coeffi¬ 
cients  relating  a  pixel's  noise  to  its  24  nearest  neighbors 
being  nonzero  (Figure  5)  , 
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Figure  5.  First  and  Second  Nearest  Neighbor 


In  Figure  5,  rk  ^  is  the  coefficient  relating  noise  components 
in  pixel  k  to  those  of  pixel  %.  Of  the  25  nonzero  values  for 
any  row  or  column,  only  six  of  them  are  distinct;  one  value 
of  1  for  the  correlation  of  a  pixel  with  itself,  four  values 
of  .3679  for  pixels  one  pixel  away,  four  values  of  .2431  for 
pixels  /2  pixels  away,  four  values  of  .1353  for  pixels  two 
pixels  away,  eight  values  of  .1069  for  pixels  /5  pixels  away, 
and  four  values  of  .0591  for  pixels  /8  pixels  away.  The 
matrix  of  correlation  coefficients  becomes 
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Using  the  simplified  noise  covariance  matrix  of  Harnly  and 
Jensen,  (6:18),  the  process  for  completing  the  generation  of 
the  noise  covariance  matrix  is  to  premultiply  the  correlation 
coefficient  matrix  by  the  variance  of  the  background  noise. 


Once  the  correlation  matrix,  R  is  known,  specific  spa¬ 
tially  correlated  realizations  of  the  64  x  1  noise  vector, 
v(t^)  can  be  computed  by  using  a  Cholesky  square  root  decom¬ 
position  of  R  and  a  zero  mean  unit-variance  independent 
Gaussian  noise  generator  (see  Appendix  E) .  The  noise  vector 
realization  will  be  added  to  the  8x8  input  array.  This 
noise  vector  is  created  by 

,v ( t^ )  =  /R  v'  (ti)  (3-32) 

where 

v' (t^)  =  64  dimensional  white  Gaussian  vector;  each 
component  zero  mean,  unit  variance,  and  in¬ 
dependent  of  other  components  (see  Appendix  E) 

?R  =  Cholesky  square  root 

The  Cholesky  square  root  decomposition  (Chapter  7  reference 
13)  produces  a  lower  triangular  matrix  such  that 

?£  /R—  =  R  (3-33) 
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By  using  this  formulation  of  the  noise  vector.  Equation 
(3-32),  the  covariance  of  v(t^)  is  shown  to  be  equal  to  the 
correlation  matrix. 

E[v(t.  )vT(t .)  ]  =  E[/|v'  (t.)v,T(t.) 

1  J  ±  J 

=  E[V' (ti)v,T(tj) ]  5lT 

=  Sr  i  ?rt  5.  . 

- -  ij 

/ 

=  R  (3-34) 


Summary  of  Truth  Model 

This  chapter  presented  the  mathematical  models  used  to 
account  for  target  dynamics,  atmospheric  jitter  effects,  and 
the  background  and  FLIR  noises.  These  models  described  the 
translational  position  changes  of  the  target's  intensity  dis¬ 
tribution  resulting  from  any  of  these  reasons.  The  mathema¬ 
tical  models  were  then  summarized  by  transforming  them  into 
state  space  notation.  Equation  (3-25)  showed  how  the  truth 
model  state  space  states  are  propagated  forward  in  time  while 
still  accounting  for  the  stochastic  nature  of  the  problem. 
Information  on  the  computer  creation  of  the  measurement  data 
for  the  computer  simulation  was  also  provided  along  with  the 
model  for  spatial  correlations  of  background  noise. 
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IV.  Extended  Kalman  Filter 


Introduction 

The  previous  chapter  presented  the  truth  model  which  is 
the  best  representation  of  the  environment  from  which  the 
measurement  vector  .z(t^)  is  generated.  An  extended  Kalman 
filter  is  used  to  process  these  measurement  vectors  to  pro¬ 
vide  an  estimate  of  the  true  target  centroid  location  for 
the  next  sample  time  (for  use  by  the  controller)  as  well  as 
of  current  state  vectors  needed  for  centering  in  the  inten¬ 
sity  image  processing.  A  four-state  extended  Kalman  filter 
model  which  accounts  for  time-corrupted  dynamics  and  the 
bandwidth  characteristics  of  the  atmospheric  jitter,  as  de¬ 
veloped  by  Mercier,  was  used  for  this  research.  A  wide 
variety  of  targets  may  be  represented  with  this  model  by 
choosing  appropriate  noise  strengths  and  correlation  times. 


Target  Dynamics  and  Atmospheric  Jitter  Filter  Models 

The  target  dynamics  and  the  atmospheric  jitter  are 
modelled  as  stationary,  first  order  Gauss-Markov  processes. 
These  processes  are  generated  as  the  outputs  of  first-order 
lags  driven  by  white  Gaussian  noise  (5:19) .  The  differential 
equation  which  describes  any  of  these  modelled  states  is 


where 


xf(t)  =  ~  xf (t)  +  wf (t) 


(4-1) 
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(4-2) 


E  [wf  (t)  ]  =  0 

E[wf (t)w.(t+x) ]  =  2of  5 (x)  (4-3) 

Tf 

and 

=  correlation  time 

wf(t)  =  white  Gaussian  noise  process 

of  »  root  mean  squared  value  of  x 

The  extended  Kalman  filter  uses  this  stationary,  first  order, 
Gauss-Markov  process  to  model  the  target  dynamics  and  the 
atmospheric  jitter  in  the  x  and  y  directions  of  the  PLIR 
image  plane.  The  stochastic  differential  equation  model  for 
each  state  upon  which  the  filter  is  based  is  Equation  (4-1) . 
This  is  the  same  structure  as  Equations  (3-2)  and  (3-3) 
except  that  a  reduced-order  model  is  used  for  atmospherics. 

State  Space  Model 

The  four  filter  states  will  now  be  put  into  state  space 

notation.  The  state  vector  differential  equation  which  de- 

T 

scribes  the  four  filter  states  (x,  x  y,  y  )  ,  is 

a  a  a  a 

^(t)  =  Fj  Xj  (t)  +  wf(t)  (4-4) 

where 

Fj  =  filter  plant  matrix 
x.j(t)  =  filter  state  vector 
w^(t)  =  input  white  Gaussian  noise  vector 
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Using  the  information  from  the  previous  section  Equation 
(4-4)  becomes 


0  0 


x^ft)  = 


+  w^t) 


0  0 


where 


0  0 


T,_  =  correlation  time  assumed  for  target  dynamics 
df 

T  _  =  correlation  time  assumed  for  atmospheric  jitter 
ar 


E  (iLf  (t)  )  =  0 


E(wf  (t)w^  (t+t)  )  =  Qf6(x) 


Sf  " 


=  assumed  target  dynamics  noise  variance 


=  assumed  atmospheric  jitter  noise  variance 


The  extended  Kalman  filter  must  propagate  its  state 
estimate  vector  and  conditional  covariance  matrix  from  one 
sample  time  to  the  next.  The  extended  Kalman  filter  model 
state  equations  are  linear  in  this  application,  and  so  the 
standard  Kalman  filter  propagation  equations  can  be  used  to 
propagate  the  filter  states  between  sample  times.  These 
standard  propagation  equations  are 

*(tI+i)  =  if(ti+i'ti)  £(ti+)  (4-9) 


£(tI+ 1>  =  if  (ti+i'ti)  i^i+)  if  (ti+i'ti)  + 


/1+1if(ti+1^>  ^  if(ti+1,r)  dx 


fci 


(4-10) 


where 


£f(ti+1#ti)  =  filter  state  transition  matrix 

P(ti+)  =  conditional  state  covariance  matrix  from 
measurement  update  equation  at  time  t^ 

P(t“+^)  =  conditional  state  covariance  matrix  pro¬ 


pagated  from  time  t^  to  t^+^ 


Qf  =  noise  covariance  kernel  descriptor  given 
in  Equation  (4-8) 

The  filter's  state  transition  matrix  $^(t^+1,ti)  must  satisfy 
the  differential  equation 

if(t,t.)  =  Sfiftt,^)  (4-id 

over  the  interval  subject  to  the  initial  condition 


— =  —  (identify  matrix)  (4-12) 

The  time-invariant  plant  matrix,  F^,  and  fixed  sample  period 
result  in  a  constant  state  transition  matrix,  (t^.^,  tA) , 
for  any  given  sample  period: 


#f(ti+i.ti)  =  e^i+l-V  =  eIAfc 


e-Afc/Tdf 


0 

0 

0 


e-^af 

0 

0 


0 

0 

e-At/Tdf 


0 
0 
0 

-At/T 


af 


(4-13) 


where  At  is  the  sample  period.  At  =  t^+^-t^. 


The  solution  to  the  integral  term  of  Equation  (4-10)  becomes 


— D= 


adf (l-e”2At/Tdf)  0 


a2fa_e-2At/Taf)  0 


a2  (l-e“2At/Tdf) 
at 


n  1 1  _-2At/T  ... 

0  a^fl-e  af) 


(4-14) 


The  matrix  of  Equation  (4-14)  is  constant  for  a  given  sampling 

time  At,  due  to  system  time  invariance  and  noise  stationarity . 
2  2 

The  values  of  and  a  ^  of  the  Qf  matrix  of  Equation  (4-8) 


45 


are  determined  during  an  off-line  tuning  process.  This  off¬ 
line  tuning  produces  the  values  which  optimize  tracking  per¬ 
formance. 

In  summary,  the  estimated  state  vector  and  the  condi¬ 
tional  covariance  matrix  propagation  equations  are 


re-At/Tdf 


x(ti+r)  = 


o 

0 

0 


e-At/Tdf 


0 

0 

0 


0 

rAt/T 

0 
0 


af 


e~At/Taf 

0 

0 


0 

-At/T 


df 


0 

0 

e-At/Tdf 

0 


0 
0 
0 

-At/T 


af 


£(t.+) 


(4-15) 
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0 
0 

-At/T 


af 


P(t.  +  ) 


e'At/Tdf 


0 

0 

0 


e-At/Taf 
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2  m _  -2At/T, 

°df  (1_e  df) 


o2f(l-e-24t/Taf>  0 


_2  / 1  -2At/T,£v 

0  cr^u-e  df) 


a2£a-e-2«/Taf» 


(4-16) 


Equations  (4-15)  and  (4-16)  propagate  the  four  filter  states 
between  sample  times.  Once  the  predicted  states  for  the 
next  sample  time  are  calculated,  the  estimates  of  the  tar¬ 
get's  position  xd(t~+1)  and  ^d^i+l^'  are  fed  to  a  feedbaclc 
controller.  The  controller  positions  the  FLIR  to  be  centered 
on  this  estimated  target  location  by  that  next  sample  time 
ti+l*  In  v;i-ew  Equation  (3-16),  the  filter  expects  the 
apparent  target's  intensity  function  to  be  offset  from  the 
center  of  that  field-of-view  by  the  estimated  jitter  effects 
also  calculated  by  Equation  (4-15) . 


Measurement  Update  Equations 

Equation  (3-30)  can  be  rewritten  in  the  general  form 


zki(ti)  ■  +  vki{ti)  (4*17) 

where  h.  . (x (t . ) , t . )  is  a  nonlinear  function  of  the  states 

KX»  X  1 

at  sample  time  t^.  Equation  (4-17)  represents  each  pixel's 
information  as  a  nonlinear  function  of  the  states  plus  addi¬ 
tive  corrupting  noise,  Vj^tt^).  Once  the  Kalman  filter's 


f 
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state  vector  and  covariance  matrix  are  propagated  to  the 
next  sample  time,  the  filter  uses  this  propagated  state  esti¬ 
mate  along  with  the  information  contained  within  each  of  the 
64  pixels  represented  by  Equation  (4-17),  to  compute  an  up¬ 
dated  estimate  of  the  underlying  target  centroid  location. 

To  process  the  information  contained  within  the  measurement 
vector,  the  extended  Kalman  filter  was  used,  since  it  can 
handle  the  nonlinearities  of  the  measurement  equation. 
Equation  (4-17) .  Moreover,  it  is  less  computationally  bur¬ 
densome  than  other  nonlinear  filters,  and  computation  time 
is  of  extreme  importance  in  this  application. 

In  the  process  of  generating  the  appropriate  filter 
gain,  an  extended  Kalman  filter  linearizes  the  measurement 
equation  about  the  most  recent  estimate  of  the  states,  x(tr)  . 
'The  inverse  covariance  form  (13:257)  of  the  extended  Kalman 
filter  measurement  update  equations  is  used  to  eliminate  the 
heed  for  the  inversion  of  a  64  x  64  matrix  for  every  update 
(5:26).  This  64  x  64  inversion  would  be  required  during  the 
calculation  of  the  Kalman  filter  gain  in  the  usual  form  of 
the  update  Equations  (13:  233),  i,e.,  K=P~HT  (HP~HT  +  R) 
because  of  the  64  scalar  measurements.  The  extended  Kalman 
filter  update  equations  in  inverse  covariance  form  are 


P_1(ti+)  =P"1(ti“)  +  HT(ti)R”1  (ti)H(t±)  (4-18) 


E(ti+)  =  [P“1(ti+)]"1 


(4-19) 
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wtrm’ 


(4-20) 


£0^)  =  £(ti+)HT(ti)R"1(ti) 

^(t±+)  =  x(ti“)  +K(ti)[z(tjL)  -htxtt/")  ,tA)  ]  (4-21) 

where 

H  (t . )  =  ^  ^  =  linearized  function  of  intensity 

1  3x  measurements  (see  Appendix  F) 

(4-22) 

P(t.~)  =  propagated  conditional  covariance  matrix  before 
1  measurement  incorporation  at  time  t. 

x(t.“)  =  propagated  state  estimate  vector  of  filter 
1  before  measurement  incorporation  at  time  t^ 

K(t^)  =  Kalman  filter  gain 

h(x(t._),t.)  =  nonlinear  function  of  intensity  measure- 
x  1  ments  at  time  t.,  as  function  of  filter 
state  estimates  x(t^'“) 

=  actual  realization  of  measurement  sector 

This  formulation  of  the  extended  Kalman  filter  update 
equations  requires  only  two  4x4  inverses,  of  the  P(t^-) 
matrix  and  the  P-1  (t^+)  •  an<*  therefore  eases  the  computa¬ 
tional  burden  of  the  algorithm  substantially.  The  inversion 
of  the  R(t^)  matrix  is  accomplished  only  once,  potentially 
offline,  and  the  result  of  that  inversion,  R-^ (t^) ,  is  stored 
for  use  in  Equation  (4-18) .  Moreover,  in  the  spatially  un¬ 
correlated  pixel  noise  case,  R-^  (t^)  is  of  the  form  (1/R)I,. 

In  this  research,  the  nonlinear  h  function  and  its 
derivatives  with  respect  to  the  filter  states  are  derived 
using  the  FFT,  phase  shifting  and  smoothing  techniques  dis- 
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cussed  in  Chapter  II.  The  structure  of  this  algorithm  is 
shown  in  Figure  1. 

In  summary,  the  inverse  covariance  form  of  the  extended 
Kalman  filter  measurement  update  equations  were  used  in  this 
research  to  ease  the  computational  burden.  Equations  (4-18) 
through  (4-21)  present  the  equations  implemented  in  the 
computer  simulation. 


Summary  of  the  Extended  Kalman  Filter  Equations 


This  section  summarizes  the  propagation  and  update 
equations  used  for  this  research 


Propagation  Equations 


x(t.-)  =  lf  x(t.^) 


(4-23 


-  if  £(ti-t>i?  +  Qd 


(4-24 


where 


e-4t/Td£ 


If' 


a-4t/Taf  0 


e-lt/Tdf  0 


e'^af 


(4-25 


and  where 


.*  y 


~D  “ 


^fa-e"2it/Tdf)  0 


2  n  — 2At/T  ... 
cf^£  (1-e  af) 


2 

Tdf' 


a^U-e 


-2At/T 


df) 


0  c2£U-e-2it'Taf> 


(4-26) 


Estimate  of  states  as  modified  by  controller  action: 


x' <t.~)  = 


0 

0 


(4-27) 


Update  Equations 


P"1  (ti+)  =  P_1  (ti~)  +  HT  (ti)  R"1^)  H(t.) 


£(ti+)  =  [P"1(ti+)]-1 


K(ti)  =  P(ti+)  HT(ti)  R1(ti) 


x(ti+)  =  x’  (tjL”)  +  K^)  [z.(ti)  -  h(x  ‘  (t^“) ,  t^)  ] 


(4-28) 

(4-29) 

(4-30) 

(4-31) 
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V,  Correlator-Kalman  Filter  Tracker 


Introduction 

The  Kalman  filter  operating  upon  the  raw  digitized 
image  has  been  shown  to  perform  well  in  comparison  to  standard 
correlation  trackers  (5; 6).  In  these  studies,  the  filter  had 
valid  apriori  information  about  the  analytic  form  of  the  in¬ 
tensity  function,  h(x(t^),t^).  The  performance  of  the  Kalman 
filter  operating  upon  the  raw  digitized  data  versus  a  cor¬ 
relator  should  be  less  beneficial  to  the  filter  without  such 
apriori  intensity  function  information.  It  is  also  computa¬ 
tionally  easier  to  implement  a  correlation  algorithm  than  to 
implement  a  high  measurement  dimension  extended  Kalman  filter. 

This  chapter  presents  an  alternate  tracking  concept 
which  is  a  hybrid  in  that  it  uses  both  correlation  and  Kalman 
filtering.  A  correlator  is  first  used  to  generate  position 
estimates  for  the  target  within  a  noise-corrupted  frame  of 
data.  The  correlator  uses  the  estimated  intensity  function, 
h(x (t^) , t^) ,  generated  via  the  data  processing  algorithm  of 
Chapter  II,  as  its  template.  This  is  different  from  current 
algorithms  which  just  use  previous  raw  data  as  a  template. 

This  modification  to  the  standard  correlation  algorithm  is 
expected  to  enhance  performance.  The  standard  correlation 
algorithm  does  not  exploit  any  knowledge  of  the  target's 
dynamics  or  any  knowledge  of  those  disturbances  which  could 
cause  apparent  translational  intensity  function  offsets. 
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The  improved  correlation  algorithm  developed  for  this  research 
positions  the  template  with  the  benefit  of  x(t^+^)  which  does 
exploit  such  knowledge.  To  account  for  correlation  errors, 
the  position  estimates  of  the  correlator  are  processed  by  a 
linear  Kalman  filter  to  produce  better  position  estimates. 

Section  one  of  this  chapter  provides  information  on  how 
the  correlation  algorithm  was  implemented.  The  next  section 
statistically  characterizes  the  errors  in  the  position  esti¬ 
mates  of  this  implementation  of  the  correlation  algorithm. 

The  last  section  of  this  chapter  develops  the  linear  Kalman 
filter  which  processes  the  correlator's  position  estimates. 


Correlator  Implementation 

This  section  presents  the  implementation  of  the  correla¬ 
tion  algorithm  used  to  provide  position  "measurements"  to  the 
Kalman  filter.  The  correlation  routine,  written  for  this 
research,  computes  the  cross  correlation  of  a  template,  the 
result  of  the  data  processing  of  Chapter  II,  and  the  noise- 
corrupted  FLIR  data  frame.  The  FFT  is  used  to  compute  the 
cross  correlation  (8:557)  as  shown  below. 


F[c[(x,y)  ] 
F[ i.(x, y)  ] 

Fta(x,y)  *  £(x,y)  ] 


G(fx,fy) 

(5-1) 

L(fX,fy) 

(5-2) 

G  (f  ,f  )  *L*(f  ,f  ) 

- —  y '  —  y' 

(5-3) 

/ 
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where 


g(x,y)  *  i(x,y)  =  cross  correlation  of  the  two-dimen¬ 
sional  spatial  sequences  g(x,y)  and 
A(x,y) 

* 

L  (f  , f  )  =  complex  conjugate  of  the  Fourier 
y  transform  of  the  sequence  Jt(x,y) 

Taking  the  inverse  transform  of  Equation  (5-3)  yields  the 

cross  correlation  of  the  two-dimensional  sequences  £(x,y) 

and  JL  (x,  y) : 

R(x,y)  =  g.(x,y)  *  A(x,y)  =  F-1  [G  (f  ,  f  )  »L  (f  ,  f  )  ]  (5-4) 

xy  Ay 

where  R(x,y)  is  the  result  of  the  correlation. 

To  show  how  the  position  estimates  are  obtained  from 

this  algorithm,  a  simple  two-dimensional  Gaussian  array  is 

examined.  Let  _fc(x,y)  represent  the  two-dimensional  array 

which  contains  the  template,  which  for  this  example  is  a 

perfect  target  replica  and  is  a  centered  Gaussian  function 
o 

with  a=2  pixels  within  the  two-dimensional  array.  Figure 
6  is  an  example  of  a  24  x  24  pixel  array  with  the  template's 
intensity  being  represented  by  a  gray-scale  (see  Subroutine 
Display  Appendix  H) .  Each  dash  around  the  perimeter  of  the 
field-of-view  represents  one  pixel  in  Figure  6.  The  use  of 
the  24  x  24  pixel  array  allows  the  8x8  pixel  tracking 
window  to  be  padded  by  8  rows  and  8  columns  of  zeros  or  data 
for  the  computation  of  FFTs . 

Similarly,  let  £(x,y)  represent  the  two-dimensional 
array  of  noise-corrupted  data,  which  contains  the  target's 


A  * 


5 


,  offset  intensity  function.  Figure  7  is  an  example  of  a  24  x 

24  pixel  array  which  is  slightly  noise  corrupted.  The  tar¬ 
get’s  single  Gaussian  intensity  pattern  is  offset  by  one  pixel 
in  the  horizontal  direction  from  the  center  of  the  24  x  24 
pixel  array.  A  gray  scale  representation  of  the  result  of 
the  cross  correlation,  R(x,y)  is  shown  in  Figure  8. 

The  expected  result  of  correlating  a  Gaussian  with  a 
Gaussian  is  a  Gaussian  whose  location  within  the  resulting 
matrix,  R(x,y),  indicates  the  relative  pixel  position  offset 
between  the  template  and  the  target.  The  two-dimensional 
array  R(x,y)  is  symmetric  in  the  vertical  direction  which 
indicates  no  offset  there.  In  the  horizontal  direction  a 
one  pixel  circular  shift  to  the  left  would  restore  symmetry. 
That  is,  if  all  of  the  columns  of  Figure  8  were  shifted  one 
column  to  the  left  and  the  leftmost  column  shifted  into  the 
rightmost  column,  horizontal  symmetry  would  be  restored. 
Obviously  the  information  that  the  target  is  one  pixel  hori¬ 
zontally  offset  from  the  center  of  the  two-dimensional  noise- 
corrupted  data  array  is  contained  in  the  result  of  the  cross 
correlation,  R(x,y).  The  fact  that  the  two-dimensional 
Fourier  transform  assumes  that  this  two-dimensional  sequence 
is  one  period  in  each  direction  of  an  infinitely  periodic 
sequence  explains  the  use  of  the  circular  shift.  Inherent 
in  the  use  of  the  FFT  to  derive  R(x,y)  is  the  assumption 
that  the  next  pixel  to  the  left  of  any  pixel  in  column  one 
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Figure  8.  Result  of  Cross  Correlation  of  Data  and  Template 
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is  equal  to  the  corresponding  pixel  in  column  24.  Using  this 
fact  to  change  the  representation  of  the  result  of  the  cross 
correlation,  the  magnitudes  of  the  pixels  intensities  of 
quadrant  one  are  switched  with  those  of  quadrant  four.  Simi¬ 
larly,  quadrant  two’s  pixels  are  swapped  with  quadrant  three's. 
The  result  of  the  quadrant  swapping  as  a  representation  of 
the  correlation  output  matrix,  R(x,y),  where  the  relative 
position  offset  between  the  template  and  the  target  is  now 
represented  by  an  offset  of  the  resulting  Gaussian's  maximum 
from  the  center  of  the  sequence.  Figure  9  provides  a  gray 
scale  representation  of  the  result  of  the  cross  correlation, 
R(x,y),  once  the  quadrants  are  swapped.  Figure  8  corresponds 
to  the  convention  of  defining  the  origin  for  the  correlation 
as  corresponding  to  the  origins  of  the  original  two  figures. 
This  means  that  the  two  figures'  origins  are  initially  super¬ 
imposed  to  create  the  first  correlation  data  point.  The 
correlation  information  corresponding  to  the  four  pixels 
which  surround  the  center  of  the  fields-of-view  of  the  ori¬ 
ginal  figures  is  contained  in  the  four  corner  pixels  in  the 
correlation  matrix.  The  convention  on  the  definition  of 
origins  attained  by  "swapping"  is  where  the  center  of  the 
24  x  24  correlation  matrix  corresponds  to  the  correlations 
of  the  centers  of  the  two  original  figures. 

By  comparing  the  original  Gaussian  template  and  the 
noise-corrupted  data  with  the  output  of  the  cross  correlation 
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given  in  Figure  9,  a  spreading  of  the  function  is  observed. 
The  noise  has  also  induced  some  asymmetries  in  the  sequence 
R(x, y) .  The  magnitudes  of  the  components  of  the  sequence 
R(x,y) ,  are  a  measure  of  the  degree  of  resemblence  between 
the  template  and  the  data  sequences.  If  the  magnitude  of  an 
element  of  the  sequence  R(x,y)  is  less  then  some  preset  frac¬ 
tion  of  the  maximum  value  of  any  element  in  that  sequence, 
then  that  element  is  considered  to  contain  poor  correlation 
information  and  can  be  set  equal  to  zero.  The  result  is  that 
such  elements  will  have  no  effect  on  the  computed  offset  be¬ 
tween  the  template  and  the  target.  Thresholding  is  often 
done  to  suppress  lower  peaks  in  the  result  of  the  correla¬ 
tions.  Figure  10  shows  the  result  of  setting  this  fraction 
to  .5  for  the  sequence  of  Figure  9. 

Once  the  threshold  has  smoothed  the  result  of  the  cross 
correlation,  a  peak  detector  is  usually  used  to  find  where 
the  maximum  resemblance  of  the  two  sequences  occurs.  Instead 
of  a  peak  detector,  a  centroid  summation  was  used  to  find 
the  center  of  mass  of  the  two-dimensional  cross  correlation 
sequence,  R(x,y).  The  center  of  mass  of  the  thresholded 
correlation  is  assumed  to  be  good  indication  of  the  peak 
location.  In  either  direction,  horizontal  or  vertical,  the 
centroid  summation  is  defined  by 


(5-5) 


N 

I  i *Amp . 


£  Amp. 
i=l  1 


where  i  is  the  horizontal  or  vertical  coordinate  of  a  given 
pixel,  and  Ampi  is  the  amplitude  value  for  that  pixel,  and  N 
is  the  total  number  of  pixels  in  the  array.  For  the  horizon¬ 
tal  direction,  the  horizontal  coordinate's  of  all  the  ele¬ 
ments  of  R(x,y)  would  be  multiplied  by  the  amplitudes  of  the 
respective  elements  and  the  products  would  be  summed.  The 
resulting  summation,  the  numerator  of  Equation  (5-5) ,  would 
then  be  divided  by  the  sum  of  all  the  amplitudes.  Equation 
(5-5)  is  implemented  in  both  directions  providing  a  position 
estimate  of  the  offset  of  the  target  from  the  center  of  the 
data  frame. 

Correlator  Error  Statistics 

This  section  statistically  characterizes  the  errors  in 
the  position  estimates  of  the  correlation  algorithm  presented 
in  the  previous  section.  These  position  estimates  will  be 
the  measurements  provided  to  a  Kalman  filter  which  will  gen¬ 
erate  a  better  estimate  of  the  target's  position.  The  two- 
dimensional  discrete-time  measurement  vector,  z(t. ),  is  a 

JL 

linear  combination  of  the  variables  of  interest,  the  target 
intensity  function's  true  position,  but  corrupted  by  an  un¬ 
certain  measurement  disturbance  yc(t^)  also  of  dimension  two. 
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(5-6) 


where 


:<V 


the  estimated  x  and  y  coordinates  by 
the  correlation/centroid  algorithm 


1100 

0011 


the  linear  combination  of  state  vari¬ 
ables  which  contribute  to  the  respec¬ 
tive  measurement  elements 


the  state  vector 


v  (t.)  =  additive  noise  corruption,  assumed  to 
be  white,  and  Gaussian,  with  statis¬ 
tics  to  be  determined 


The  additive  noise  corruption  vector,  V£ (t^)  of  the  measure¬ 
ment  model  of  Equation  (5-6),  must  account  for  the  statistical 
effects  of  the  errors  in  the  correlator/ centroid  position 
estimates . 

Software  was  developed  to  test  the  correlator/ centroid 
algorithm's  position  estimates  repeatedly,  and  histograms  of 
the  resulting  error  were  generated.  Figure  11  contains 
examples  of  the  error  histograms  produced  for  a  target  off¬ 
set  by  .3  pixels  from  the  template  in  the  horizontal  direc¬ 
tion  only.  The  errors  could  be  modelled  as  Gaussian  random 
variables  with  appropriate  means  and  variances.  The  values 
needed  to  describe  the  error's  mean  or  variance  would  be  a 
function  of  the  template-target  separation,  the  threshold 
level  used  by  the  centroid  and  even  the  variance  of  the 
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Figure  11a.  Histogram  of  Error  in  Horizontal  Position  Estimate 


Figure  lib.  Histogram  of  Error  in  Vertical  Position  Estimate 


background  noise.  For  this  research  a  value  for  the  threshold 
of  .5  and  a  target-template  separation  of  .15  pixels  were 
chosen  to  establish  error  statistics.  Other  values  which  were 
considered  for  the  threshold  were  .1,  .2,  .25,  .75,  and  even 
.90.  The  value  of  .5  was  chosen  since  it  produced  acceptable 
results  for  the  signal-to-noise  ratios  of  interest.  Error 
data  for  other  choices  of  thresholds  will  be  provided  in  the 
next  chapter.  At  the  sampling  rate  of  30Hz,  for  the  trajec¬ 
tories  being  simulated  in  this  research,  a  propagation  error 
of  not  more  than  .15  pixels  was  expected.  Background  noise 
was  allowed  to  vary  over  the  full  range  of  signal-to-noise 

ratios  of  10  to  20.  The  result  was  a  mean  error  of  approxi- 

2 

mately  -.080  pixels  and  a  variance  of  .00363  pixels  in  the 

horizontal  direction.  The  error  in  the  vertical  direction 

2 

has  a  mean  of  .116  pixel  with  a  variance  of  .00598  pixels  . 
These  values  can  now  be  used  to  describe  statistically  the 
errors  in  the  correlator/centroid  position  estimates.  This 
characterization  is  needed  for  the  measurement  model  of 
Equation  (5-6) . 

The  difference  in  the  vertical  and  horizontal  errors 
above  is  because  of  the  locations  of  the  three  Gaussians  of 
the  target  intensity  profile.  If  the  center  of  the  8x8 
pixel  tracking  window  is  located  at  coordinates  (0,0)  then 
the  Gaussians  would  be  located  at  (0,-2.667),  (-2,1.333), 

and  (+2,1.333).  This  placement  forces  the  center  of  mass  of 
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the  intensity  profile  to  coordinates  (0,0),  but  the  resulting 
intensity  pattern  is  not  radially  symmetric  about  (0,0). 

Kalman  Filter  Tracker 

The  mathematical  models  for  target  dynamics  and  atmos¬ 
pheric  jitter  which  were  developed  in  Chapter  IV  are  also 
used  in  the  Kalman  filter  which  processes  the  correlator's 
intensity  function  position  estimates.  These  models  accounted 
for  time  correlated  dynamics  and  the  bandwidth  characteristics 
of  the  atmospheric  jitter.  The  standard  Kalman  filter  propa¬ 
gation  equations,  as  given  by  Equations  (4-23)  through  (4-26), 
are  also  used  by  this  Kalman  filter.  The  difference  in  the 
Kalman  filter  used  in  this  alternative  tracker  is  that  it  pro¬ 
cesses  a  two-dimensional  measurement  vector  of  intensity 
function  position  estimates,  whereas  the  previous  filter  pro¬ 
cessed  64  intensity  measurements.  The  measurement  equation 
for  the  new  filter  was  given  in  Equation  (5-6) .  The  low 
dimensionality  of  the  measurement  vector  eliminates  the  need 
and  in  fact  the  applicability  of  the  inverse  covariance  form 
of  the  Kalman  filter  update  equations.  The  correlator/cen¬ 
troid  intensity  function's  position  estimates,  which  consti¬ 
tute  the  two-dimensional  measurement  vector,  are  linear  com¬ 
binations  of  the  Kalman  filter  states,  as  seen  in  Equation 
(5-6; .  The  linear  filter  measurement  update  equations  used 
are  given  in  Equations  (5-7)  through  (5-9) . 


(5-7) 


K(t^)  =  P(ti")HT(ti)  [H(ti)P(ti")HT(ti)  +R(ti)]_1 

xtt^)  =  x(t..~)  +K(ti)[z(ti)  -  a(t±)x(t”)  ]  (5-8) 

£(ti+)  =  P(ti~)  -  K(ti)H(ti)P(ti')  ( 5— S ) 

The  values  of  P(t“)  and  x(ti~)  are  obtained  from  the  propaga¬ 
tion  equations  and  the  measurement  uncertainty  covariance 
matrix,  R (t^) ,  is  assumed  diagonal  with  the  values  determined 
as  explained  in  the  previous  section.  The  diagonal  nature  of 
R(t^)  assumes  that  the  correlation  position  estimate  uncer¬ 
tainty  in  one  direction  is  independent  of  the  uncertainty  in 
the  position  estimate  in  the  other  direction.  Since  the  errors 
are  a  function  of  SNR  and  threshold  level,  this  independence 
assumption  is  subject  to  later  revision.  Direct  computation 
of  the  off-diagonal  elements  of  R(ti)  could  verify  or  invali¬ 
date  the  assumption  of  independence  of  the  errors  in  the  po¬ 
sition  estimates. 

In  summary,  this  chapter  developed  an  alternate  tracking 
algorithm  which  uses  a  correlation  algorithm  to  provide  inten¬ 
sity  function  position  estimates  as  a  measurement  vector  to  a 
Kalman  filter.  This  algorithm  is  different  from  standard  cor¬ 
relation  trackers  in  that  it  uses  the  model  of  target  dynamics 
from  the  Kalman  filter  to  position  the  template,  thresholding 
to  remove  false  peaks,  on-line  derived  template  shape,  and 
correlator  position  estimate  enhancement  via  a  Kalman  filter. 
The  Kalman  filter  uses  internal  models  to  separate  what 
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translational  motion  of  the  intensity  function  resulted  from 
target  dynamics  and  what  resulted  from  atmospheric  jitter. 
The  results  from  this  tracking  algorithm  and  the  algorithm 
of  Figure  1  are  given  in  the  next  chapter. 


t 
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VI.  Performance  Analysis 

Introduction 

This  chapter  presents  the  results  from  the  testing  of 
the  tracking  algorithms  developed  in  the  previous  chapters. 
The  first  section  of  this  chapter  derives  figures  of  merit 
which  provide  a  means  to  evaluate  the  accuracy  in  the  deri¬ 
vation  of  the  intensity  functions  contained  in  Chapter  II. 
This  accuracy  is  ultimately  only  important  in  its  impact  on 
tracking  ability.  The  sensitivity  of  tracking  performance  to 
the  accuracy  of  the  derived  intensity  functions  will  be  shown 
as  well.  The  tracking  ability  performance  criteria  are  pre¬ 
sented  in  the  second  section  of  this  chapter.  The  third 
section  discusses  those  variable  parameters  which  the  com¬ 
puter  simulation  allows  the  input  to  change  for  sensitivity 
analysis.  The  errors  in  the  derivation  of  the  intensity 
functions  and  the  tracking  errors  are  plotted  as  mean  +  one 
sigma  (standard  deviation)  errors  in  the  next  section.  The 
final  section  of  this  chapter  condenses  the  results  of  this 
research  into  tables.  These  tables  cross  reference  varia¬ 
tions  in  the  parameters  which  control  the  Kalman  filter  and 
the  pattern  recognition  process  to  tracking  and  intensity 
function  derivation  errors.  Important  trends  and  sensitivi¬ 
ties  are  revealed  in  these  tables.  The  statistical  informa¬ 
tion  of  this  chapter  was  generated  using  Monte  Carlo  tech¬ 
niques  (13:329),  see  Appendix  G. 
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Derivation  of  Intensity  Functions 

The  accuracy  of  the  on-line  derivation  of  the  intensity 
functions  (see  Chapter  II)  is  actually  of  importance  here 
only  since  it  affects  tracking  ability.  To  make  some  reason¬ 
able  judgement  on  how  good  a  representation  of  h(x(t^),t^) 
and  H(jc(t^),t^)  has  been  derived,  a  simple  figure  of  merit  is 
developed  in  this  section. 

At  each  frame  the  truth  model  provides  information  on 
the  location  of  the  centroid  of  the  true  target's  intensity 
functions.  The  data  processing  of  Chapter  II  derives  on-line 
estimates  of  these  intensity  functions.  At  each  frame  the 
difference  between  the  true  target  intensity  functions: 

1.  the  nonlinear  intensity  function,  hfc 

2.  the  derivative  of  ht  with  respect  to  a  change  in 
either  horizontal  state,  xd  or  xg, 

3.  the  derivative  of  h^  with  respect  to  a  change  in 
either  vertical  state,  yd  or  ya,  Hyt 

A  A  A 

and  the  estimated  intensity  functions,  h,  H  ,  and  H  ,  respec- 

~x  ~y 

tively,  are  computed  pixel  by  pixel. 

As  the  simulation  progresses,  at  each  frame  i  two  8x8 
arrays  are  maintained  to  compute  the  accuracy  of  .  One 

A 

array  is  the  difference  between  ht^  and  h^.  The  square  of 
this  error  is  also  maintained:  (hti  -  h^) 2  where  i  represents 
which  frame  is  being  processed. 

Many  time  histories  are  processed  in  the  Monte  Carlo 
study  (as  described  in  detail  in  Appendix  G) ,  and  for  each  j 
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frame  the  mean  and  variance  of  the  error  associated  with  each 
pixel  for  that  frame  is  desired.  The  mean  error  for  a  given 
pixel,  j  =  1,  .  .  . ,  64,  for  a  given  time  frame,  i  =  1,  .  .  ., 
20,  over  all  the  time  histories  run,  k  =  1,..,N,  would  be 
given  by 

,  N  „  ,  N 

(ht  "  hiik}  n  E  e—  (6-D 

ijk  1JK  ,  k=l 


1  N 

Mo  =  £  Z 
eij  N  k=l 


Sijk 


where 


M  =  mean  error  for  pixel  j  at  the  time  frame  i, 
eij  averaged  over  the  N  simulations. 

N  =  number  of  simulations 

k  =  index  of  Monte  Carlo  simulation  runs,  going 
from  1  to  N 

h.  =  value  for  frame  i,  pixel  j,  of  the  truth  model 
rijk  intensity  for  the  kth  simulation 


h.  ..  =  estimate  value  for  frame  i,  pixel  j,  of  the 
1-'  intensity  as  derived  in  the  kth  simulation 

The  variance  of  the  error  for  frame  i,  pixel  j,  would  be 

given  by 


2  1  N 

Te  =  N  E 
eij  W  k=l 


1  N 
“  N  E 


k=l 

N 


(eijk  Meij) 

(eL*  '  2Me:  eijk  +  » 


1J 


2M 


=  I  v  e- 

N  "=1  eijk  N  ~ij 

,  N 


N 

Z 

k=l 


'ij 

1  N 

eijk  +  N  e.j 


„  „  e2  -  2M 
N  k=1  ijk  e.j 


Z 

k= 
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M  +  “*N »M  2 
6ij  N  6ij 
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Note  that  in  general  1/(N-1)  could  be  used  in  place  of  1/N 
in  the  first  line  of  this  equality  to  reduce  the  bias  in  this 
estimate,  but  with  N=20  this  substitution  would  not  make  much 
difference.  Equations  (6-1)  and  (6-2)  show  that  two  8x8 
arrays  must  be  maintained  for  each  frame  for  each  intensity 
function  of  interest.  In  this  research,  the  errors  in  fi,  H^, 
and  5^  are  of  interest,  and  so  six  arrays  were  maintained  for 
each  frame. 

For  this  research,  new  data  frames  are  generated  at  a 
30  Hz  rate.  If,  as  in  Appendix  G,  twenty  frames  constitute 
a  single  tracking  time  history,  two-thirds  of  a  second  simu¬ 
lation  time,  then  six  8  x  8  x  20  arrays  are  needed.  At  the 
end  of  the  simulation  these  six  8  x  8  x  20  arrays  contain 
the  information  needed  to  compute  the  ensemble  average  and 
variance  of  errors  for  any  pixel,  j,  for  any  frame,  i.  with 
the  assumption  that  every  pixel  is  as  important  as  every 
other  pixel,  the  spatial  average  of  the  mean  pixel  errors 
can  be  computed.  Equation  (6-3)  shows  how  to  compute  the 
spatial  average  of  the  mean  pixel  errors  for  any  frame,  i. 
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j=l 


(6-3) 


where 


=  the  spatial  average  (over  all  the  pixels)  of 
ei  the  mean  pixel  errors  at  the  ith  frame 

M  =  mean  (i.e.  ensemble  average  over  all  the  simu- 
eij  lations)  at  the  jth  pixel  for  the  ith  frame 
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Note  that  if  each  pixel  were  not  considered  equally  important, 
a  weighted  average  could  replace  (6-3) .  Similarly,  the  spa¬ 
tial  average  of  the  variance  for  any  frame  i  could  be  com¬ 
puted  by 


(6-4) 


where 


2 

oe  =  the  spatial  average  of  all  the  ensemble  pixel 
l  error  variances,  over  all  the  pixels  at  the 
ith  frame 


=  variance  of  the  error  of  the  jth  pixel  at  the 
ith  frame 


The  benchmark  for  intensity  function  derivation  errors  is 

A  A 

found  by  providing  perfect  information  of  xpeak  and  Tpeak  to 
the  pattern  recognition  algorithm.  This  benchmark  is  dis¬ 
played  as  the  twelfth  entry  of  Table  I  in  the  next  chapter. 


Tracking  Ability 

Appendix  G  describes  the  Monte  Carlo  study  used  to  gen¬ 
erate  the  statistics  of  this  section.  The  quantities  of 
interest,  with  respect  to  tracking,  are  the  errors  in  the 
estimated  value  of  x^(t7) ,  yd(t~),  xd(t^),  and  It 

is  also  important  to  estimate  the  location  of  the  intensity 
distribution's  centroid.  The  accuracy  of  the  estimation  of 
the  centroid's  location,  xpeak(ti+p  and  ypeak(ti+l)  of 
Equation  (3-30a) ,  will  affect  the  accuracy  of  the  estimated 
intensity  functions  and  in-turn  affect  the  tracking  accuracy. 
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By  comparing  the  errors  before  and  after  incorporation  of  a 
measurement,  information  about  how  the  estimates  were  improved 
from  the  information  of  the  data  frame  is  obtained. 

Figure  G-l  shows  how  the  error  samples  are  generated. 
Similar  to  what  was  done  in  the  previous  section,  error  sta¬ 
tistics  are  calculated  as  the  mean  error  for  the  variable  of 
interest  at  frame  i. 


(6-5) 


where 

E  =  mean  error  (i.e.  ensemble  average  error  over 
xd^  all  simulations)  in  x  dynamics  for  frame  i 

=  truth  model,  x-dynamics  value  at  frame  i  for 
dik  simulation  k 

A 

x.  =  filter  estimated  x-dynamic  value  at  frame  i 
aik  for  simulation  k 

and  the  variance  of  the  errors , 


1  N 
=  -  l 


<3. 

l 


2  —  2 
„  -  e  -  E 

N  k=l  xdik  x 


(6-6) 


di 


The  1/ (N-l)  substitution  discussed  in  the  previous  section 
is  also  applicable  here.  The  generalization  of  equations 
(6-5)  and  (6-6)  to  compute  the  errors  in  yd,  Xpealc/  or  ypea^ 
is  clear.  The  benchmark  for  tracking  performance  is  set  by 
providing  perfect  h,  H^,  and  H^.  This  is  the  case  studied 
by  Capt  Mercier  (5).  Table  II  of  that  research  (5:57)  lists 
a  mean  error  of  .2  pixels  for  a  SNR  =  20.  That  value  can  be 
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compared  to  entry  one  of  Table  I,  contained  in  the  next 
section,  to  show  only  a  slight  degradation  in  tracking 
ability  when  deriving  the  intensity  profiles. 

Variation  of  Parameters 

The  computer  simulation  was  developed  to  allow  certain 
filter  and  data  processing  parameters  to  be  varied.  This 
was  accomplished  in  order  to  test  the  tracking  algorithms  of 
Figure  1  and  Chapter  V  with  different  design  characteristics 
and  in  different  tracking  environments. 

The  first  parameter  to  be  varied  is  the  spread  para- 
2 

meter,  ag  ,  of  the  target's  Gaussian  intensity  profiles. 

This  value  is  used  by  the  truth  model  to  generate  the  three- 
Gaussian  hot  spot  data  each  having  the  given  spread,  (see 
Equation (3-35) ) .  The  standard  value  for  this  parameter,  for 
this  research,  was  2.0  but  runs  were  also  made  with  values 
of  1.0  and  3.0. 

The  next  parameter  is  the  number  of  zeros  to  pad  around 
the  data.  Padding  of  the  finite  data  array  with  zeros  is  a 
common  engineering  practice  to  allow  manipulation  of  the 
transformed  data  array  to  provide  results  compatible  with 
the  limited  f ield-of-view.  In  Chapter  II  it  was  shown  that 
DFTs  assume  that  the  finite  data  array  is  one  period  of  an 
infinitely  periodic  two-dimensional  sequence.  Any  manipula¬ 
tion  of  the  transformed  data  array  uses  this  assumption  of 
infinite  periodicity.  Padding  with  8  rows  and  8  columns  of 


zeros  creates  a  24  x  24  pixel  array  which  is  assumed  to  be 
one  period  by  the  DFT.  The  padding  insures  that  the  infinite 
periodicity  assumption  will  not  affect  results  within  the 
8x8  tracking  window.  If  the  spread  parameter,  ag  ,  is 
chosen  such  that  the  target  intensity  height  is  approaching 
zero  near  the  edge  of  the  8x8  tracking  window,  then  padding 
with  zeros  will  not  adversely  affect  the  results  of  the  data 
processing  of  Figure  1.  However,  if  ag  is  such  that  signi¬ 
ficant  intensity  magnitudes  exist  outside  the  8x8  array, 
to  pad  with  zeros  arbitrarily  would  induce  an  artificial  edge 
in  the  intensity  function.  These  edges  in  the  two-dimensional 
spatial  intensity  array  will  cause  increased  magnitudes  of 
the  high  frequency  components  in  the  transform  domain.  For 
this  application,  a  full  frame  of  FLIR  data  actually  consists 
of  300  x  400  pixels,  so,  when  necessary,  the  8x8  array 
could  be  padded  with  the  noise-corrupted  data  instead  of  zeros . 
When  this  parameter  is  set  to  eight,  then  the  8x8  tracking 
window  is  surrounded  by  8  rows  and  8  columns  of  zeros  to  fill 
up  a  24  x  24  complex  data  array.  The  standard  value,  for 
this  research,  for  this  parameter  is  zero  causing  the  8x8 
tracking  window  to  be  padded  by  noise-corrupted  data. 

The  next  two  parameters  are  the  number  of  frames  per 
simulation  and  the  number  of  simulations  per  Monte  Carlo 
study.  Appendix  G  explains  these  parameters  in  terms  of  the 
Monte  Carlo  study.  Both  of  these  parameters  were  set  to  20. 


78 


I 


Early  in  the  testing  of  the  tracking  algorithms,  fifty  frames 
were  used  to  insure  steady-state  error  was  reached  in  a  20 
frame  time  history.  The  consistency  of  the  computed  statis¬ 
tics,  between  twenty  and  fifty  frame  time  histories  motivated 
the  use  of  the  twenty  frame  simulation. 

Alpha,  the  relative  weighting  parameter  for  the  smooth¬ 
ing  process,  is  the  next  parameter  which  may  be  varied. 
Chapter  II  explains  how  alpha  affects  the  smoothing  process, 
and  Equation  (2-9)  explicitly  displays  its  role  in  the  algo¬ 
rithm.  For  this  research,  the  standard  value  for  alpha  is 
.1,  but  a  value  for  alpha  of  1.  and  .2  are  also  investigated. 
Using  alpha  equal  to  1.  results  in  no  smoothing  and  thus 
provides  a  benchmark  to  demonstrate  what  benefit  smoothing 
provides . 

The  number  of  high  frequency  components  of  the  FFT  of 
the  image  to  zero  out  is  the  next  parameter.  This  provides 
the  ability  to  investigate  how  spatial  filtering  within  the 
intensity  function  derivation  portion  of  Figure  1  could  en¬ 
hance  or  corrupt  tracking  performance.  Spatial  filtering 
could  easily  be  accomplished  within  the  optical  implementa¬ 
tion  discussed  in  the  next  chapter  if  it  is  shown  beneficial 
here.  For  standard  runs,  no  frequencies  were  zeroed.  Runs 
zeroing  the  two  and  four  highest  spatial  frequency  components 
were  accomplished  also. 

2 

The  input  background  noise  variance,  ,  is  next. 


*  * 
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This  value  also  includes  FLIR  noise  contributions.  Since 
the  maximum  values  of  the  three  Gaussian  intensities  were 
twenty,  the  signal  to  noise  ratio,  SNR,  is  equal  to 

SNR  =  atak-signal  value - 20  ( 

rms  background  noise  a ^  ' 

SNR  values  of  20  (standard)  to  10  were  considered  as  repre¬ 
sentative  of  realistic  tracking  scenarios. 

The  next  two  parameters  are  filter  parameters  which 
are  varied  to  tune  the  filter  for  optimal  tracking.  The 
variance  of  the  dynamic  discrete  time  noise  driving  the  tar¬ 
get  states  of  the  filter  (see  Equation  4-8)  is  a  measure  of 
the  urcertainty  of  the  filters  dynamics  model.  The  variance 
of  the  atmospherics  for  the  filter  (see  Equation  4-8)  is 
similarly  explained.  Without  any  in-depth  tuning  to  optimize 
tracking  performance,  values  of  .1  and  .001  pixels2  were  used 
respectively.  The  disparity  in  these  values  indicates  the 
filter's  atmospheric  model  is  expected  to  be  3  very  good 
representation  of  the  atmospheric  disturbance.  The  approxi¬ 
mation  of  ignoring  the  high  frequency  pole  made  by  the  Kalman 
filter  is  not  expected  to  affect  the  filter's  atmospheric 
representation  to  any  substantial  degree  based  on  performance 
analyses  of  previous  filters  that  incorporated  this  same 
approximation  (4? 5; 6).  These  values  did  establish  acceptable 
tuning.  Chapter  VIII,  Recommendations,  discusses  issues  of 
tuning  which  should  eventually  be  considered. 


The  next  parameter  is  the  RMS  atmospheric  jitter  in¬ 


duced  in  pixels  (see  Equation  3-8).  A  value  of  .1  pixel 
jitter  was  used  as  a  standard. 

The  software  also  allows  the  option  of  which  domain 
the  exponential  smoothing  algorithm  is  to  be  accomplished. 
Appendix  C  discusses  this  algorithm  in  both  domains.  Fre¬ 
quency  domain  smoothing  is  used  for  standard  runs  in  this 
research. 

Plotting  Results 

Plots  of  the  errors  in  the  algorithm's  representation 
of  the  intensity  functions  and  in  tracking  are  presented  in 
this  section.  For  each  setting  of  the  parameter  values  of 
the  previous  section,  15  plots  were  generated.  All  of  the 
plots  are  mean  errors  +  one  standard  deviation. 

The  first  four  plots  are  of  the  errors  in  the  estimates 
of  the  x  and  y  location  of  the  target.  Errors  at  a  minus 
time  are  the  errors  before  incorporation  of  a  measurement 
while  errors  at  a  plus  time  are  after  the  measurement  has 
been  processed  at  that  time.  Plots  five  and  six,  of  any  15 
plot  sequence,  are  the  errors  in  x  or  y  dynamic  estimates  at 
a  minus  time  followed  immediately  by  the  errors  at  the  plus 
time.  Observing  the  way  these  plots  are  driven  toward 
smaller  errors  at  plus  times  relative  to  minus  time  estimates 
clearly  shows  how  much  information  is  obtained  from  that 
frame.  Recall  that  the  dynamics  state  estimates  are  needed 
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for  the  controller  and  errors  in  these  estimates  are  the  fun¬ 
damental  indicator  of  tracking  performance,  while  centroid 
estimates  are  required  for  centering  images  in  the  upper 
processing  path  of  Figure  1.  Plots  seven  through  twelve  are 
the  corresponding  error  plots  for  the  centroid  estimates. 
Plots  13  through  15  are  plots  of  the  intensity  function  spa¬ 
tial  average  derivation  errors  of  Equations  (6-3)  and  (6-4) . 

For  each  of  the  fifteen  plots,  a  legend  is  provided  to 

show  how  the  variable  parameters  of  the  previous  section  were 

2 

set  for  a  given  plot.  This  legend  consists  of  the  COV  =  a  , 

g 

NZ  =  number  of  zero  to  pad,  rows  and  columns  used,  ALP  =  a 
used  for  smoothing,  NF  =  number  of  frequencies  to  zero,  VAB  = 
variance  of  the  background  and  FLIR  noise,  and  SDF  =  sigma  of 
the  dynamics  assumed  by  the  filter.  The  following  15  plots. 
Figures  12  through  27,  are  provided  here  as  an  example  of  the 
sequence  and  results,  for  the  case  of  COV=2,  NZ=0,  ALP=.l, 
NF=0,  VAB=1.,  SDF= . 1 .  Appendix  I  contains  plots  for  four 
other  settings  of  the  variable  parameters. 

Figures  12  and  13  show  the  x  and  y  errors  at  a  minus 
time.  The  plots  indicate  that  the  steady  state  error  for 
either  variable  is  reached  by  the  fourth  frame.  The  sigma 
is  approximately  constant  after  that  point  also.  The  same 
trend  is  also  seen  in  the  x  and  y  plus  time  errors,  Figures 
14  and  15,  except  that  the  errors  are  smaller  as  expected 
after  incorporation  of  a  measurement.  Plots  16  and  17  show 
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how  much  information  is  obtained,  and  the  corresponding  re¬ 
duction  in  errors,  by  incorporating  a  measurement.  For  these 
plots  the  errors  at  a  minus  time  are  plotted  just  prior  to 
their  corresponding  plus  time  errors.  It  is  clear  from  these 
plots  that  consistently  good  information  is  being  obtained 
from  the  noise-corrupted  measurements  under  the  standard 
parameter  settings.  Figures  18  through  23  provide  similar 
plots  for  the  centroid  location  estimates.  The  same  trend 
is  evident  for  these  plots  except  that  the  errors  are  smaller 
This  result  was  expected  in  that  it  is  easier  to  find  cen¬ 
troids  than  it  is  to  separate  dynamic  and  atmospheric  contri¬ 
butions  to  intensity  function  movement.  Plots  24  through  27 
present  the  intensity  function  derivation  errors.  These 
plots  clearly  show  that  mean  errors  are  reduced  very  quickly 
to  very  small  values?  the  standard  deviations  are  also  re¬ 
duced  accordingly.  The  inadequacy  of  the  dynamics  models 
used  by  the  Kalman  filter  resulted  in  the  propagated  state 
estimates  being  consistently  biased  in  the  same  direction  in 
these  plots.  Since  it  was  the  main  goal  of  this  research  to 
develop  and  implement  tracking  algorithms  which  derive  the 
target  intensity  functions  in  an  on-line  manner,  the  adequacy 
of  the  dynamics  model  was  not  emphasized.  This  problem  is 
readily  rectified  by  a  dynamics  model  more  representative  of 
the  target  characteristics  than  a  simple  first  order  Gauss- 
Markov  position  process  in  each  direction.  Moreover,  the 
propagation  errors  also  provide  a  means  to  determine  how 
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Figure  12.  X  Minus  Errors 
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Figure  16.  X  Position  Error:  xd(tA  )  and  xd(ti  )  Errors 
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Figure  17.  Y  Position  Error:  Yj(tj  )  and  yd(tt  )  Errors 
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Figure  20.  X  Centroid  Pius  Error 
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MEAN  ERROR  +/-  SIGMA  OF  ESTIMATED  INTENSITY  PROFILE 


much  information  is  being  obtained  during  measurement  updates 
when  these  on-line  derived  intensity  functions  are  being  used. 


Table  Summary 

This  section  will  summarize  the  tracking  and  intensity 
function  estimation  abilities  in  tabular  form.  The  plots  of 
the  previous  section  show  errors  plotted  as  a  function  of 
frame  number.  The  results  displayed  in  Tables  I-III  are  the 
corresponding  time  average  of  those  errors  and  their  standard 
deviations  over  frames  10  to  20.  The  plots  of  the  previous 
section  implied  that  the  steady  state  error  region  is  reached 
by  this  time,  and  thus  time  averaging  is  confined  to  this  in¬ 
terval.  The  special  comment  column  is  to  designate  pertinent 
information  not  indicated  by  the  header  (e.g.  smooth  in  space, 
run  made  on  Cyber  175/Modcomp  minicomputer) . 

The  first  result  from  Table  I  is  that  exponential 
smoothing  produces  similar  results  in  either  the  frequency  or 
space  domain.  Entries  two,  four,  seven,  and  nine  of  Table  I 
are  compared  to  the  entries  directly  above  them  to  show  this 
result.  This  was  expected  from  the  derivation  of  exponential 
smoothing  contained  in  Chapter  II.  There  are  additional 
errors  induced  when  smoothing  in  space  is  used  for  this  im¬ 
plementation.  Recall  Figure  1  of  Chapter  I.  The  upper  path 
of  that  figure  is  where  smoothing  is  accomplished.  Before 
smoothing  can  be  accomplished,  the  intensity  function  must 
be  centered  within  the  data  frame.  The  centering  is 


Table  I.  Tracking  and  Intensity  Function  Derivation  Errors 
for  Modcomp  (minicomputer)  (32  bit) 


(COV,NZ,ALP,NF,VAB,SDT) 


Table  II.  Tracking  and  Intensity  Function  Derivation  Errors 
for  Cyber 
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Table  III.  Tracking  Errors  for  Correia tor-Kalman  Filter 
for  Modcomp 
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accomplished  in  the  frequency  domain,  (see  Chapter  II) .  If 
smoothing  is  to  be  accomplished  in  the  space  domain,  the 
centered  data  must  then  be  inverse  transformed.  The  smoothed 
data  would  then  be  transformed  again  for  evaluation  at  the 
propagated  state  estimate  and  differentiation  which  are  accom¬ 
plished  in  the  frequency  domain.  The  additional  transforms 
required  for  smoothing  in  space  induce  additional  numerical 
errors.  The  very  slight  difference  between  the  results  shows 
that  transforming  errors  are  practically  negligible.  Smooth¬ 
ing  in  space  may  be  used  for  the  optical  implementations 
which  will  be  presented  in  the  next  chapter. 

Table  I  contains  the  results  obtained  on  the  Mod comp 
while  Table  II  contains  the  results  obtained  on  the  Cyber 
175.  The  Modcomp  minicomputer  uses  32  bit  arithmetic  com¬ 
pared  to  the  Cyber's  60  bits.  The  first  entry  from  either 
table  is  the  standard  run  as  defined  by  the  parameter  settings 
contained  in  the  header.  The  comparable  results  obtained  in 
tracking  errors  for  example  in  *err(+)'  -.03399  +  .10559  com¬ 
pared  to  -.06009  +  .10155,  reveal  an  insensitivity  of  the 
developed  algorithms  to  word  length.  Similarly,  comparable 
results  were  obtained  for  filtering  out  high  frequencies, 
padding  with  zeros,  and  even  the  maximum  noise  corruption 
case.  This  result  is  important  since  it  implies  an  insensi¬ 
tivity  to  h,  H  ,  and  H  errors.  The  insensitivity  also  im- 
x  y 

plies  that  h  and  the  Kalman  filter  gain,  which  uses  H  and 
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Zeroing  out  the  two  highest  spatial  frequencies  is  seen 

to  increase  accuracy  on  all  columns  in  Table  I.  This  was 

2 

expected  for  the  spread  variable,  ,  equal  to  two.  (This 
will  be  explored  further  when  broad  and  sharp  intensity 
shapes  are  discussed  below.)  The  values  for  the  three 
Gaussian  dispersions,  whose  peaks  are  separated  by  approxi¬ 
mately  four  pixels,  results  in  an  intensity  function  which 
contains  relatively  low  spatial  frequencies.  The  high  spa¬ 
tial  frequency  components  are  because  of  noise  and  the  zero¬ 
ing  of  them  decreases  the  intensity  function  derivation  errors 
and  thereby  also  increases  tracking' accuracy.  Tables  I  and 
II  show  that  zeroing  out  the  four  highest  spatial  frequencies 
didn't  result  in  a  significant  change  in  the  errors  from  the 
two  frequency  cut  off  case. 

Padding  with  zeros  slightly  increases  errors  over  the 
pad  with  data  case  for  the  standard  parameter  settings.  This 
shows  that  the  true  intensity  profile  is  not  equal  to  zero 
outside  the  8x8  tracking  window  for  the  spread  parameter 
equal  to  2.  To  pad  with  zeros  introduces  artificial  high 
frequencies  as  explained  earlier,  and  zeroing  out  these  arti¬ 
ficial  high  frequencies  improved  tracking  by  18  percent  for 
x(-)  (see  Table  II  entries  5  and  6),  where  only  a  2  percent 
improvement  was  achieved  when  padding  with  data  (see  Table  II 


rr.  **.’**•' 


entries  1  and  2).  For  this  application  the  luxury  of  padding 
with  data  exists  because  a  frame  of  FLIR  data  actually  con¬ 
sists  of  300  x  400  pixels.  Padding  with  zeros  will  be 
shown  to  enhance  tracking  for  sharp  intensity  profile  cases 
below.  For  sharp  target  intensity  profiles  the  true  target 
intensity  is  approximately  zero  outside  the  8x8,  so  only 
noise  contributes  significant  amplitudes  there.  Zeroing 
that  noise  enhances  tracking. 

The  eighth  entry  of  Table  I  provides  the  tracking  and 
intensity  function  derivation  errors  when  SNR  is  at  its 
minimum  level  used  in  this  research,  i.e.,  SNR  =  10.  The 
increased  noise  did  degrade  tracking,  by  about  40%  for  x(-) 
and  by  about  300%  for  x(+),  as  was  expected  since  it  becomes 
harder  to  find  the  target  within  the  data  frame.  The  deriva¬ 
tion  errors  for  h,  for  the  maximum  noise  corruption  case  of 
Table  II  entry  seven,  is  seen  not  to  change  significantly 
from  the  case  of  high  standard  SNR  as  in  entry  one.  This 
result  shows  that  the  pattern  recognition  algorithm  is  not 
very  sensitive  to  the  increase  in  noise. 

Entries  four,  seven,  and  nine  of  Table  I  present  in¬ 
formation  on  the  combination  of  smoothing  in  space  with  cut¬ 
ting  the  two  highest  spatial  frequencies,  padding  with  zeros 
or  in  combination  with  maximum  noise  corruption.  Smoothing 
in  space  for  all  of  these  cases  produced  results  compatible 
with  smoothing  in  the  frequency  domain.  This  is  seen  when 


each  of  the  smoothing  in  space  entries  are  compared  with  the 
entries  directly  above  them  in  Table  I. 

Entry  ten  of  Table  I  and  entry  sixteen  of  Table  II 
result  from  using  an  alpha  equal  to  .2  instead  of  the  standard 
.1.  Tracking  ability  is  degraded  from  that  of  entry  one  even 
though  errors  in  the  derivation  of  h  and  Hy  are  not  increased. 
However,  errors  in  deriving  do  increase.  Since  setting 
alpha  equal  to  .2  simulates  a  finite  memory  filter  of  about 
five  samples  long,  this  result  shows  that  a  higher  magnitude 
steady  state  error  level  is  reached  for  Hx  because  less 
smoothing  is  accomplished.  But,  if  the  target  shape  does 
change,  this  larger  alpha  would  respond  faster  than  the 
smaller  standard  alpha.  Recall  that  the  placement  of  the 
three  Gaussians,  as  explained  in  Chapter  V,  is  the  reason 
that  Hx  is  harder  to  derive  than  Hy.  The  no  smoothing  entry 
of  Table  II  shows  that  no  smoothing  increases  tracking  errors, 
x(+),  by  almost  400%.  Cutting  four  frequencies  reduces 
tracking  errors  to  about  300%  of  the  smoothed  error. 

The  sharply  peaked  intensity  entries  of  Table  II  show 
that  tracking  and  intensity  function  derivation  errors  in¬ 
crease  over  the  corresponding  broad  intensity  cases.  The 
increased  errors  are  because  the  target's  intensity  functions 
are  not  as  easily  distinguished  from  the  noise  as  when  the 
target  has  a  broad  intensity  profile.  Padding  with  zeros 
for  the  sharply  peaked  intensity  entries  of  Table  II  decreases 
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tracking  errors  by  10*.  Recall  that  padding  with  zeros  in- 

O 

creased  errors  for  the  a *  =  2  case.  The  placement  of  the 
Gaussians  with  small  dispersions  result  in  the  target  only 
contributing  insignificant  magnitudes  outside  the  8x8 
tracking  window.  The  padding  with  zeros  only  essentially 
zeros  only  noise  for  these  cases  resulting  in  smaller  errors. 

Last  it  should  be  noted  that  the  Correlator-Kalman 
filter  algorithm  of  Chapter  V  outperformed  the  Figure  1  algo¬ 
rithm  consistently.  The  first  entry  of  Table  III  shows  that 
for  the  standard  case  that  the  correlator-Kalman  filter  algo¬ 
rithm  has  smaller  tracking  errors,  25*  less  for  x(+),  and 

centroid  location  errors,  32*  less  for  x  .  (+) .  The  sensi- 

peak 

tivity  of  this  algorithm  to  the  variation  of  parameters  is 
similar  to  the  results  discussed  above.  Padding  with  zeros 
increased  errors,  by  35*  for  x(-),  while  cutting  high  fre¬ 
quencies  did  not  have  much  effect  on  the  standard  case. 

Maximum  noise  increases  errors,  by  20*  for  x(+),  but  cutting 
high  frequencies  from  this  minimum  SNR  case  regains  9*  of 
that  lossed  accuracy.  Setting  alpha  equal  to  .2,  entry  8 
Table  III,  is  found  not  to  induce  additional  tracking  errors, 
unlike  the  effect  it  had  on  the  algorithm  of  Figure  1.  Since 
this  algorithm  does  not  use  Hx,  this  result  was  expected. 

Also  it  must  be  re-emphasized  that  this  is  not  a  standard 
correlator.  This  correlation  algorithm  uses  dynamic  modelling 
information  from  the  Kalman  filter  to  position  the  template. 
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Thresholding  is  used  to  avoid  false  peaks  along  with  increas¬ 
ing  the  accuracy  of  the  centroid  calculation.  Thresholds  of 
.1  and  .9  were  tested  for  no  noise  and  for  SNRs  of  10  and  20. 
The  mean  error,  for  SNR  =  10,  decreased  from  .08211  pixels 
for  no  threshold  to  -.03870  pixels  for  a  threshold  of  .2. 
Although  these  mean  errors  are  both  close  to  zero  and  their 
difference  may  be  considered  insignificant,  the  standard  de¬ 
viation  made  a  significant  decrease  from  .01357  to  .00531 
pixels.  At  a  threshold  of  .5  the  mean  error  had  increased 
to  .080  pixels  and  the  standard  deviation  had  dropped  to 
.00363  pixels.  The  histogram  of  the  errors  also  approximated 
the  Gaussian  form  more  than  the  .2  or  .3  threshold  cases. 

This  Gaussian  form  along  with  the  tighter  standard  deviation 
led  to  the  choice  of  using  a  threshold  of  .5  for  this  research 
The  third  enhancement  of  the  correlator  is  that  the  template 
is  being  derived  on-line  via  the  data  processing  of  Chapter 
II.  Finally,  the  enhanced  correlator  position  estimates  are 
processed  by  a  Kalman  filter.  Capt  Mercier's  research  (5) 
tested  a  standard  correlation  tracker  against  an  extended 
Kalman  filter  algorithm  which  utilized  apriori  knowledge  of 
the  intensity  functions.  In  Table  II  of  that  study  (5:57) 
the  mean  track  error,  for  a  standard  correlation  tracker,  is 
listed  at  .5  pixels  with  a  la  error  of  1.5  pixels  for  a 
SNR  =  20.  The  correlation  algorithm  developed  for  this  re¬ 
search  resulted  in  a  mean  error  of  -.02543  pixels  with  a  la 


error  of  .13437  pixels.  This  algorithm  shows  considerable 
potential  as  a  next  generation  tracker. 

Summary 

This  chapter  presented  the  results  from  the  testing  of 
the  tracking  algorithms  developed  in  the  previous  chapters. 
The  errors  are  plotted  as  mean  +  one  sigma  (standard  devia¬ 
tion)  errors  and  tables  of  errors  are  generated  for  variable 
settings  for  Kalman  filter  and  pattern  recognition  para¬ 
meters.  Extremely  good  tracking  is  established  for  both  al¬ 
gorithms  with  the  Correlator-Kalman  filter  algorithm  showing 
considerable  potential. 

The  better  performance  of  the  alternate  tracking  algo¬ 
rithm  can  be  explained  as  follows.  Previous  comparisons 
(4 ; 5)  of  Kalman  filter  trackers  to  correlators  used  apriori 
knowledge  of  the  intensity  functions.  As  less  knowledge  is 
available,  the  processing  of  the  64  intensity  measurements 
via  a  Kalman  filter  is  less  accurate,  with  less  relative 
benefit  over  performance  attainable  with  a  correlator.  More 
over,  as  discussed  in  the  previous  section,  this  alternate 
tracker  embodies  much  more  than  a  simple  correlation  tracker 
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Introduction 


I 

A 


The  computational  burden  of  either  of  the  tracking 
algorithms  developed  in  the  previous  chapters  is  heavy.  For 
that  reason,  a  serious  attempt  needs  to  be  made  to  develop 
hybrid,  optical  and  digital,  systems  for  implementation  of 
those  algorithms.  Optical  processing  offers  the  potential 
of  parallel  processing  of  two-dimensional  data.  The  two- 
dimensional  Fourier  transforms,  correlations,  and  matrix 
operations  which  are  contained  in  the  tracking  algorithms 
could  potentially  be  accomplished  at  extremely  high  data 
rates  using  optical  techniques.  A  digital  computer  could 
be  programmed  to  perform  necessary  data  analysis  or  control 
the  optical  computations.  Such  a  hybrid  system  would  com¬ 
bine  the  speed  and  parallel  processing  advantages  of  optical 
processors  with  the  flexibility  of  digital  computers. 

This  chapter  begins  with  a  section  which  describes  the 
basic  principles  of  optical  processing  and  introduces  the 
concept  of  hybrid  processing.  The  next  section  adapts  the 
hybrid  processor  to  the  tracking  algorithms  developed  pre¬ 
viously.  Potential  modifications  to  the  tracking  algorithms 
are  proposed  to  make  optical  implementation  easier.  The 
purpose  of  this  section  is  to  present  design  possibilities 
and  to  provide  an  extremely  rough,  preliminary  analysis  of 
the  potential  for  fruition  of  optical  implementations  for 
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each  proposal. 


Background 

This  section  presents  the  basic  principles  of  optical 
processing  systems.  These  systems  utilize  optical  interfer¬ 
ence  to  process  two-dimensional  incoming  signals.  The  most 
common  system  used  for  optical  processing  is  shown  in  Figure 
27.  Light  from  a  laser  point  source  S  is  collimated  by  lens 
Lc*  The  collimated  beam  of  quasi-monochromatic  coherent 
laser  light  is  used  to  illuminate  the  object  in  plane  P-^  and 
a  diffraction  pattern  is  formed.  The  object  is  placed  one 
focal  length  in  front  of  the  transforming  lens  L-^.  The  class 
ical  Fraunhofer  diffraction  pattern  of  the  object's  space- 
varying  amplitude  transmittance,  t^x^y^,  in  plane  is 
formed  in  the  rear  focal  plane  of  lens  .  The  complex  field 
of  the  Fraunhofer  diffraction  pattern  in  plane  P2  is  mathema¬ 
tically  equal  to  the  spatial  Fourier  transform  (7; 166)  of  the 
object. 

T0  (x2,y2>  =  kj/Zt^x^y^expl  ^-{x^  +  y^)  Idx-jdyj^  (7-1) 

T0(fx'fy)  =  ki//to(xi,y1)exp[-j2lT{fxx1  +  fyy1)]dx1dy1  (7-2) 

where 

k^  =  complex  constant 

T 0  =  Fourier  transform  of  tQ(x1#y^) 

X  =  mean  wavelength  of  the  laser 
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Figure  27.  Basic  Optical  Processing  System 


The  two-dimensional  Fourier  transform  of  the  object  in  plane 
P^L  can  be  modified  by  placing  an  appropriate  optical  filter 
in  plane  P 2»  This  optical  filter  is  in  general  a  complex 
function.  The  optical  filter  diffracts  the  light  incident 
upon  plane  P^,  producing  a  second  diffraction  pattern  (14:28). 
This  diffraction  pattern  is  then  focussed  by  lens  L2  onto 
plane  P^ . 

In  a  hybrid  system,  the  input  plane  P^ ,  the  transform 
plane  P2,  and  the  output  image  plane  P^  are  all  interfaced 
via  spatial  light  modulators  to  a  digital  computer.  To  uti¬ 
lize  the  pro<" ->ssing  capability  of  the  optical  processor,  some 
means  of  inputting  information  into  and  reading  information 
out  of  these  spatial  light  modulators  at  data  rates  commen¬ 
surate  with  the  30Hz  FLIR  frame  repetition  rate  is  required. 
Spatial  light  modulators  with  cycle  frequencies  of  30Hz  are 
now  available  from  several  vendors. 

Itek  Corporation,  for  example,  produces  a  Pockels  Read¬ 
out  Optical  Modulator,  (PROM) .  The  PROM  is  a  spatial  light 
modulator  which  is  constructed  from  bismuth  silicon  oxide 
(bisox) .  Bisox  is  a  crystal  which  displays  both  wavelength 
selective  photoconductivity  and  the  linear  electro-optic 
(Pockels)  effect  (15:353).  The  complexity  of  the  total  ex¬ 
posure-readout  process  has  precluded  the  development  of  a 
complete  and  exact  theoretical  model  of  i ROM  operation  for 
predicting  system  response  for  varying  operating  conditions 
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(15:354).  However,  simplified  models  which  are  limited  to 
particular  phases  of  the  process  have  been  developed.  These 
theoretical  models  include  predictions  on  the  wavelength 
dependence  of  PROM  sensitivity,  the  relationship  between  PROM 
sensitivity  and  Bisox  crystal  thickness,  and  how  PROM  resolu¬ 
tion  can  be  increased  by  increasing  capacitance  of  the  block¬ 
ing  layer  (15:364).  The  characteristics  of  the  PROM,  or  of 
any  of  the  spatial  light  modulators  considered  in  this  analy¬ 
sis  do  not  lend  themselves  to  presentation  in  a  simple  tabular 
form  because  of  the  many  variables  involved  in  the  complex 
exposure-readout  process.  However,  certain  benchmarks  on 
commercially  available  PROMs  are  available.  A  25mm  diameter 
is  standard  with  special  order  devices  with  37mm  diameters 
also  available.  The  desired  30Hz  cycle  rate  is  standard  with 
a  dc  contrast  ratio  of  103  to  104:1.  A  spatial  frequency 
response  of  100  cycles  per  millimeter  has  been  demonstrated 
(15:364)  . 

Hughes  also  produces  a  spatial  light  modulator  which 
could  be  used  in  this  application.  The  hybrid  field-effect 
liquid  crystal  light  value  (LCLV)  is  a  high  resolution  opti- 
cal-to-optical  image  converter.  The  input  and  output  light 
beams  are  completely  separated  and  noninteracting.  The  para¬ 
meters  which  characterize  LCLV  performance  are:  sensitometry; 
modulation  transfer  function  and  resolution;  linearity; 
signal-to-noise  ratio;  response  time;  optical  flatness;  and 
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image  quality  (16:374).  The  performance  of  the  LCLV  can  only 
be  given  for  specific  values  of  the  operational  parameters 
(voltage  bias  and  frequency,  and  orientation  of  the  light 
value  with  respect  to  the  output  light  polarization  direc¬ 
tion)  .  The  performance  is  particularly  sensitive  to  the 
imaging  light  power.  The  use  of  intensifier  tubes  which  are 
capable  of  amplifying  very  low  light  level  scenes  to  the 
levels  required  by  the  light  value  is  feasible  (16:381). 

These  are  not  the  only  available  spatial  light  modula¬ 
tors.  Their  characteristics  are  however,  indicative  of  cur¬ 
rently  available  devices. 


Applying  Optical  Processing  to  Tracking 

This  section  will  apply  the  optical  processing  tech¬ 
niques  discussed  in  the  previous  section  to  the  two  tracking 
algorithms  developed  in  chapters  one  through  six.  This  sec¬ 
tion  will  also  propose  potential  modifications  to  the  track¬ 
ing  algorithms  to  facilitate  hybrid  implementations. 

The  data  processing  algorithm  of  Chapter  II  is  first 
considered  for  optical  implementation.  Both  tracking  algo¬ 
rithms  use  this  procedure  for  generating  information  about 
the  target’s  intensity  profile.  In  addition,  the  algorithm 
of  Figure  1  also  uses  the  derivative  information  developed. 
The  software  implementation  of  the  pattern  recognition  algo¬ 
rithm,  Appendix  H,  along  with  the  explanation  of  Chapter  II 
clearly  shows  how  manipulation  of  the  amplitudes  and  phases 


of  the  various  spatial  frequencies  are  used  to  derive  the 
respective  intensity  functions.  A  complex  filter  placed  in 
the  transform  plane,  P of  Figure  27  can  be  used  to  manipu¬ 
late  the  amplitudes  and  phases  of  those  spatial  frequencies 
simultaneously.  The  spatial  filtering  results  of  the  pre¬ 
vious  chapter  become  important  here.  Complex  filters  are 
usually  implemented  using  holographic  techniques  which  were 
pioneered  by  A.  B.  Vander  Lugt  of  the  University  of  Michigan's 
Radar  Laboratory.  The  best  hybrid  combination  for  generation 
of  the  flexible  complex  filter  which  would  accomplish  the 
application  of  the  appropriate  linear  phase  shift  to  center 
the  target's  intensity  function  within  the  incoming  data 
array  requires  further  research.  Once  this  filter  or  combin¬ 
ation  of  filters  is  synthesized  and  centering  of  the  data  is 
accomplished,  using  the  configuration  of  Figure  27  with  the 
complex  filter  discussed  above  in  plane  P plane  P^  would 
contain  the  space  domain  representation  of  the  centered  tar¬ 
get  intensity  function.  The  next  step  is  to  smooth  this 
image  with  the  previous  best  estimate  of  the  intensity  func¬ 
tion  stored  in  a  spatial  light  modulator.  Since  exponen¬ 
tially  smoothing  can  be  accomplished  in  either  domain,  as 
established  in  Chapter  VI,  there  is  no  constraint  on  where 
an  optical  implementation  performs  it.  Both  images  are  mul¬ 
tiplied  by  constants  (see  Chapter  II)  and  the  products  are 
added,  probably  via  simultaneous  recording  techniques,  to 
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form  the  new  best  estimate  of  the  intensity  function  which 
is  written  into  the  spatial  light  modulator  which  stores  the 
smoothed  estimate.  This  smoothed  estimate  must  also  be  pro¬ 
cessed  by  a  complex  shifting  filter  to  produce  the  intensity 
function  expected  at  the  next  sample  time,  h(x(ti),ti).  Dif¬ 
ferentiation  of  this  expected  intensity  profile  can  be  accom¬ 
plished  using  a  simple  amplitude  filter  in  the  transform 
plane,  P2  of  Figure  27.  The  derivative  of  the  object,  in 
plane  of  Figure  27  can  be  produced  by  a  filter  function 
whose  amplitude  transmittance  is 

f(kx)  =  a(kQ  -  kx)  (7-3) 

where 

a,kQ  =  constants 

The  use  of  the  bias  term  kg  eliminates  the  need  for  a  complex 
filter  (15:361)  for  differentiation.  This  amplitude  trans¬ 
mittance  filter  has  its  maximum  transmittance  at  the  center 
of  the  f ield-of-view  with  a  minimum  transmission  on  the 
edges.  This  is  analogous  to  using  the  differentiation  prop¬ 
erty  of  the  Fourier  transform  as  explained  in  Chapter  II. 
Digitally  controlled  complex  filters  could  theoretically 
accomplish  the  Chapter  II  intensity  function  derivations. 

Such  filters  are  not  easily  implemented.  The  information 
which  is  to  be  input  into  the  transform  plane  filter  spatial 
light  modulator  would  be  computed  by  a  digital  computer. 

The  computed  desired  transfer  function  is  added  to  the 
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reference  wave  within  the  digital  computer  to  simulate  the 
holographic  synthesis  technique.  The  output  of  the  digital 
computations  is  a  non-negative,  real-valued  function  which 
is  quantized.  Since  the  desired  masks  are  of  finite  extent, 
recording  of  sampled  values  of  the  computed  transfer  function 
can  completely  specify  its  Fourier  transform.  This  assumes 
that  the  well  known  sampling  theorem  (7:25)  is  satisfied. 

The  major  advantage  of  digitally  constructed  filters  is  that 
an  abstract,  mathematically  defined  filter  can  be  computed 
and  recorded.  This  eliminates  the  need  for  the  availability 
of  the  point-spread  function.  Such  complex  digital  computa¬ 
tions  could  minimize  the  advantage  of  a  hybrid  implementation 
for  a  small  tracking  window  such  as  the  8x8  pixels  used  in 
this  research.  However,  if  this  intensity  function  could  be 
shown  to  be  only  slowly  changing  a  slower  repetition  rate  on 
its  generation  could  relieve  the  digital  computational  burden 
Another  possible  application  of  optical  processing  to 
the  algorithm  of  Figure  1  would  be  in  the  implementation  of 
the  high  dimensioned  matrix  operations.  The  parallel  nature 
of  optical  processing  could  be  used  to  great  advantage  in 
accomplishing  these  functional  operations.  The  formation  of 
the  expected  intensity  function,  h(x(t^),t^),  was  discussed 
above.  Once  the  spatial  domain  representation  of  Mx(t^),t^) 
is  created,  it  is  to  be  subtracted  from  the  measurement  data 
z.(t.),  which  is  the  incoming  FLIR  data  frame.  The  intensity 
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values  from  the  FLIR  could  drive  a  spatial  light  modulator 
and  then  the  residual  could  be  formed  by  optical  subtraction 
of  jz(t^)  -  h(x(t^),t^)  (16:383).  The  subtraction  scheme  of 

Reference  16,  page  383,  uses  two  LCLVs  onto  which  z.(t^)  and 
h(x(t^,ti)  are  projected.  Potentially  the  IR  radiation  which 
was  input  to  the  FLIR  could  directly  drive  a  spatial  light 
modulator  eliminating  the  need  for  a  scanning  type  FLIR.  The 
residual  must  next  be  multiplied  by  the  gain  matrix  Kft^) . 

In  the  digital  implementation  this  matrix  is  4  x  64  with 
only  two  unique  rows.  Therefore,  two  two-dimensional  images 
which  represent  the  two  unique  rows  could  be  placed,  via 
spatial  light  modulators,  in  the  plane,  in  the  space  domain, 
which  contains  the  residual  resulting  in  two-dimensional 
images  which  could  be  sampled  and  summed  to  produce  the  two 
unique  update  components  of  Equation  (4-2) .  These  are  the 
components  which  are  added  to  our  best  estimate  of  the  state 
before  we  get  a  measurement,  x (t~) ,  to  generate  x(t^+) . 

To  enhance  the  algorithm  of  Chapter  V,  Correlator- 
Kalman  filter,  optical  correlations  can  be  used.  The  trans¬ 
form  of  the  template  will  be  encoded  into  the  filter  of  the 
transform  plane,  of  Figure  27.  This  template  could 

theoretically  be  obtained  by  accomplishing  the  data  processing 
of  Chapter  II  optically  as  discussed  above.  A  second  possi¬ 
bility  would  be  to  derive  the  template  digitally  as  was  done 
for  Chapter  V  and  then  input  the  information  into  the  spatial 
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light  modulator;  this  was  also  discussed  above.  A  third 
possibility  would  be  to  array  many  potential  templates  in 
the  filter  plane  and  accomplish  simultaneous  optical  correla¬ 
tions  with  all  of  them.  Templates  could  be  prepared  in  which 
the  presence  and  disposition  of  targets  were  varied  and  where 
the  background  characteristics  were  also  varied.  The  weight¬ 
ing  of  the  results  from  the  multiple  template  correlations 
could  be  controlled  via  adaptive  estimation  of  target  and 
environmental  characteristics .  The  input  data  image  would 
be  used  to  modulate  a  collimated  laser  beam  which  is  input 
to  a  transform  lens.  A  multiple  holographic  lens  is  used  as 
the  Fourier  transforming  lens  which  directs  the  transform  to 
each  of  the  many  matched  filters,  made  from  the  templates. 

The  position  estimates  from  any  of  these  correlation  imple¬ 
mentation  options  could  then  be  processed  by  a  Kalman  filter. 

Conclusion 

In  summary,  high  bandwidth  optical  processing  in  par¬ 
allel  with  flexible  lower  bandwidth  digital  computation 
potentially  could  relieve  the  computational  burden  of  the 
tracking  algorithms  developed  in  this  research.  The  pattern 
recognition  algorithm  given  in  Chapter  II  would  be  very  dif¬ 
ficult  to  implement  optically  because  of  the  complex  filters 
required.  For  implementation  of  the  high  dimensional  Kalman 

filter,  (see  Chapter  IV) ,  optical  processing  could  accomplish 

* 

the  computationally  intense  functional  matrix  operations. 
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Optical  implementation  of  the  correlation  algorithm  of  Chap¬ 
ter  V  has  the  most  promise.  Before  any  optical  implementation 
is  possible  though,  considerable  research  and  development  must 
be  accomplished.  The  software  developed  in  this  research 
provides  a  flexible  tool  to  determine  what  spatial  frequency 
modifications  enhance  tracking.  The  zeroing  out  of  spatial 
frequencies  within  the  software  is  just  one  example  of  how 
this  software  can  help  determine  what  filters  are  needed  for 
optical  implementation  of  the  tracking  algorithms  and  what 
characteristics  the  resulting  hybrid  optical/digital  algorithm 
would  have.  The  testing  of  the  differentiation  filter  of 
Equation  (7-3)  would  just  be  a  point  by  point  multiplication 
of  the  spatial  frequency  components  by  an  array  whose  values 
were  determined  by  that  equation.  The  software  provided  in 
Appendix  H  is  well  commented  and  easily  modified  for  this 
type  of  testing. 


Conclusions 


The  conclusions  addressed  in  this  section  are  formulated 
from  the  details  of  Chapter  VI,  Performance  Analysis.  The 
milestones  which  were  realized  in  this  research  will  be  reit¬ 
erated  here. 

As  shown  in  Table  I,  the  extended  Kalman  filter  algo¬ 
rithm  performs  very  well  even  when  the  intensity  functions 
are  not  assumed,  as  was  done  in  prior  research  efforts  (4; 5). 
These  intensity  functions  were  derived  using  the  digital 
pattern  recognition  techniques  developed  in  Chapter  II.  It 
was  shown  in  Chapter  VI  that  the  tracking  capability  was  not 
overly  sensitive  to  increased  errors  in  the  derivation  of 
the  intensity  functions.  The  filtering  of  high  spatial  fre¬ 
quencies  was  shown  to  decrease  tracking  errors  for  cases 
where  the  intensity  function  is  relatively  smooth. 

In  previous  research  efforts,  apriori  intensity  func¬ 
tion  information  was  given  to  the  extended  Kalman  filter. 

For  some  tracking  scenarios  those  extended  Kalman  filter 
trackers  were  found  to  outperform  correlation  trackers.  The 
alternate  tracking  algorithm  of  Chapter  V  was  developed  to 
combine  correlation  and  Kalman  filtering.  When  no  apriori 
information  on  the  intensity  functions  was  assumed,  the  pro¬ 
cessing  of  the  64  intensity  measurements  by  the  Kalman  filter 
should  not  outperform  correlators  so  easily.  For  this  reason 
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a  correlation  algorithm  which  makes  use  of  dynamic  model 
information  from  the  Kalman  filter,  thresholding  correlator 
outputs  to  remove  false  peaks  prior  to  centroid  calculations, 
on-line  template  derivation  and  correlator  position  estimate 
enhancement  via  a  Kalman  filter  was  derived.  Chapter  VI 

1 

I 

showed  how  this  algorithm  outperforms  the  originally  envi¬ 
sioned  algorithm  of  Figure  1. 

Finally  Chapter  VII  presented  the  optical  processing 
alternatives  for  implementing  the  tracking  algorithms.  The 
parallel  nature  of  optical  processors  provide  the  potential 
to  relieve  the  computational  burden  of  the  tracking  algo¬ 
rithms.  Prior  to  any  optical  implementation,  extended  re¬ 
search  and  development  must  be  accomplished. 

Recommendations 

Further  research  is  recommended  to  assess  the  effects 
of  incorporating  the  IR  laser's  FLIR  image  in  the  model's  j 

used.  The  use  of  optical  windows  could  potentially  notch 
the  contribution  of  the  IR  laser  out  of  the  FLIR  image.  The 
laser  will  however  induce  hot  spots  on  the  target  which  were 
not  part  of  the  original  multi-hotspot  target  model.'.  Modern 
estimation  techniques  could  also  be  used  to  determine  what 
translational  or  phase  distortion  effects  occurred  during 
propagation  of  the  laser  beam.  The  use  of  spatial  light  modu-  , 

lators  could  potentially  obtain  the  phase  distortion  induced 

j 

by  the  atmosphere.  A  compensating  imposed  distortion  on  the  : 
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outgoing  beam  could  then  maximize  target  damage. 

Another  area  of  suggested  research  is  to  attempt  opti¬ 
cal  implementation  of  those  portions  of  the  tracking  algo¬ 
rithms  discussed  in  Chapter  VII.  This  must  include  research 
into  the  feasibility  of  using  the  many  types  of  available 
spatial  light  modulators. 

The  alternative  tracking  algorithm  developed  in  Chapter 
V  needs  to  be  investigated  further.  The  optimal  thresholding 
and  peak  detection  or  centroiding  techniques  could  further 
improve  tracking  ability.  Adaptive  techniques  for  on-line 
tuning  could  also  be  investigated. 

The  algorithms  of  this  research  should  now  be  tested 
using  real  FLIR  data.  A  combined  Physics-EE  effort  at  AFIT 
would  be  required  to  assemble  the  experimental  apparatus 
needed  to  generate  the  real  FLIR  data.  As  an  alternative, 
the  data  could  potentially  be  provided  to  AFIT  by  the  Air 
Force  Weapons  Laboratory. 

The  use  of  different  dynamics  models  should  also  be 
investigated.  Adaptive  models  or  models  compatible  with 
different  tracking  scenarios  could  also  be  investigated. 

An  investigation  of  how  these  algorithms  perform  when 
tracking  a  slowly  changing  target  shape  is  also  needed.  The 
variation  of  alpha  to  enhance  tracking  of  slowly  changing 
targets  could  also  be  investigated. 

Finally,  only  empirical  convergence  of  the  algorithms 
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developed  was  shown  during  this  research.  A  theoretical 
proof  of  convergence  of  these  adaptive  estimation  techniques 
should  be  derived. 
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jendix  A 


Two-Dimensional  Finite  Discrete  Fourier  Transform 
Interpretation  and  Test 


Similar  to  the  presentation  of  J.  W.  Goodman,  reference 
7,  the  Rect  function  is  defined  to  extract  one  period  of  the 
spatial  intensity  function  as 


RectN(x,y)  =  1  0<x<N-l,  0<y£N-l 

0  otherwise 


(A-l ) 


Reproducing  Equations  (2-3)  and  (2-4)  of  Chapter  II  and  using 
the  Rect  function  to  define  the  area  sequence  g(x,y)  as  being 
zero  outside  the  interval  0<x<N-l 


"x  y 


N-l 

N-l 

Z 

Z 

x=0 

1! 

o 

g(x,y)  =  g  (x,y)  RectN(x,y) 


z,l?X  (fx  +  f  Y) 
m  x  y 


(A-2) 


Rectff  ,f  ) 
N  x  y 


N-l  N— 1 


+  j2TT  (f  x  +  f  v) 


g(x,y)  =  ~  Z  Z  G  (f  ,  f  )  e 

tr  f  =o  f  =o  y 

x  y 


(A-3) 


Recwy 

(A-4 ) 


G(f  ,f  )  is  the  Finite  Discrete  Fourier  Transform  of  g(x,y) 
x  y 

(8:118) .  The  Rect  function  of  Equation  (A-l)  is  separable 
in  the  two  independent  variables: 


Rect„(f  ,  f  )  =  Rect..  (f  )  RectM(f  ) 
N  x  y  N  x  N  y 


Equation  (A-3)  can  now  be  written  as 


(A-5) 


-  i  2tt  (f  y ) 


G  (f  / f  )  =  Z  X{fx,y)exp 
*  y=0 


RectN(fy) 


(A-6) 


where 


N_!  -J2tt  (fxx) 

X(f  ,y)  =  I  g  (x,y) exp  N  Rect  (f  ) 

x=0 


(A-7 ) 


Equation  (A-7),  X(fx,y)  corresponds  to  an  N-point  one  dimen¬ 
sional  Discrete  Fourier  Transform  for  each  value  of  the  row 
index  y  (8:118).  X(fx,y)  is  the  result  of  N  one  dimensional 
transforms,  one  for  each  row  of  g(x,y).  In  Figure  A-la,  if 
y  is  held  constant,  say  y=2,  and  a  one-dimensional  Fourier 
Transform  is  accomplished  across  that  row,  the  variation  in 
the  image  intensities  from  column  to  column  would  result  in 
the  zero  frequency  of  that  variation  to  be  located  at  x'=l 
and  y'=2,  of  Figure  A-lb,  while  the  coefficient  associated 
with  the  fundamental,  the  first  harmonic,  in  both  directions 
would  be  found  at  x'=2  and  y'=2  with  its  conjugate  located 
at  x'=N  and  y'=2. 

To  complete  the  two-dimensional  Discrete  Fourier  Trans¬ 
form  Equation  (A-6)  expresses  how  to  implement  the  N  one¬ 
dimensional  transforms  along  each  column.  To  summarize  this 
general  interpretation  of  the  two-dimensional  Discrete  Fourier 
Transform,  the  Equations  (A-6)  and  (A-7)  show  how  the  two- 
dimensional  transform  is  achieved  by  using  a  one-dimensional 
transform  on  the  rows  first  and  then  on  the  columns,  or  vice 
versa. 


Figure  A-l.  Row  Transform  Result  of  2-D  DFT 


For  the  case  where  the  image  intensity  function  is 
separable  in  the  two  independent  variables  x  and  y,  for  ex¬ 
ample,  g(x,y)  having  the  property  that 

g(x,y)  =  g1(x)g2(y)  (A-8) 

the  two-dimensional  Discrete  Fourier  Transform,  G(f  ,f  ), 

x  y 

becomes  the  product  of  the  one-dimensional  independent  trans¬ 
forms  G^f^  and  G2(fy). 

G(fx,fy)  =  G1(fx)G2(fy)  (A- 9 ) 

A  function  of  two  independent  variables  is  separable  within 
a  specific  coordinate  system  if  it  can  be  written  as  a  pro¬ 
duct  of  two  functions  of  one  independent  variable  each.  The 
two-dimensional  transform  degenerates  to  holding  the  row  index 
constant,  for  example,  and  running  the  one-dimensional  trans¬ 
form  across  the  columns,  for  anyone  of  the  N  rows.  Similarly 
the  column  index  is  held  while  the  one-dimensional  transform 
is  accomplished  across  the  rows,  and  this  is  done  for  any  of 
the  columns.  The  results  of  these  two  independent  transforms 
are  then  multiplied  together. 

This  leads  to  some  interesting  algebra  which  can  be 
exploited  to  test  the  implementation  of  the  two-dimensional 
Discrete  Fourier  Transform.  In  Figure  A-l,  T^2  2^  will  con¬ 
sist  of  the  product  of  the  two  one-dimensional  transform 
fundamental  coefficients  for  that  row  or  column.  If  x  is 
the  column  coordinate  and  y  is  the  row  coordinate,  then 
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T (  xj  is  defined  as 


(2,2) 


Is  harmonic  in  y  =y 
for  column  2 


(1,2) 


1  harmonic  in  x  =x 


(2,N) 


1st  harmonic  in  y  =y,,  fc 

for  column  N  U'N)  *  h 

j  J  [Lf 


for  row  2 


conjugate  of  V 


(1,2) 


(A-10) 


harmonic  in  x  =x 
for  row  2 


(1,2) 


(A-ll) 


T(N-1,2)  =  [conjugate  _  of  1st 


harmonic  in  y  =y 
for  column  2 


T  (N— 1  ,N)  =  Icon  jugate  of  lSt|  A 


(1,  2) 


1  harmonic  in  x  =x 
for  row  N-l 


(1 , N-l) 


(A-12) 


for  column  N 


harmonic  in  y  j  =y  ^  Ny*  |*  ]  jharmonic  in  x  |  =x 


conjugate  of  1  A  * 


for  row  N-l  J 


(1 , N-l 


(A-13) 


if  y1<N/2,  xx<N/2 

Tty^x-^  =  fliy!-] 


-1)  harmonic  in  y  =y,  ,  , 

for  column  x,  vyl  ,xl' 


(x,-l)  harmonic  in  x  =x 


1  for  row  yx  ^•Y1) 

If  the  variation  in  intensity  across  every  row  is  equal  to 
the  variation  in  intensity  for  every  other  row  and  the  vari¬ 
ation  across  all  of  the  columns  is  similarly  set  equal,  then 


the  components  of  the  transform  can  be  defined  as 


=  a  +  jb 


y(l,2) 
yU,N) 
y  (1  /  2) 
y  (1  #  N) 


=  a  +  jb 


*  =  a  -  jb 

*  =  a  -  jb 


x(l,  2)  °  °  +  Jd 

X  (1,  N— 1 )  '  C  +  Jd 
XU,2)*  "  c  '  Jd 
XU,N-1>  “  c  - 


(A-14) 


Using  Equation  (A-14)  to  solve  (A-10)  through  (A-13) 

T (2  2)  =  Y(1  2)  *x  (i  2)  =  (a+Jb)  (C+Jd)  =  (ac-bd)  +  j  (bc+ad) 

(A-15) 


T(2,N)  =  y  (1 , N)  *X*  (1 , 2)  =  (a+Jb)(c-jd)  =  (ac+bd)  -  j  (ad-bc) 

(A-16) 

T  (N-l ,  2)  =  *(1,2)  X(1,N-1)  =  (a"jb)  *  {c+jd)  =  <ac+bd>  +  J  (ad-bO 

(A-17) 

T  (N-l ,  N)  -  y(l,N)xa,N-l)  =  (a-jb)  •  (c-jd)  =  (ac-bd)  -  j  (ad+bc) 

(A-18) 

Equation  (A-15)  for  2)  tbe  conJu9ate  of  Equation  (A-18) 

for  t(n-1  n)  under  the  conditions  of  uniform  variations  on 
any  row  and  similarly  for  the  columns,  which  just  implies 
separability. 

To  implement  the  two-dimensional  Fourier  Transform, 
subroutine  Fourt  was  used,  which  is  a  common  Fortran  sub¬ 
routine  (4:).  To  test  subroutine  Fourt  to  see  where  it  places 
harmonics  with  a  two-dimensional  array,  the  results  of  Equa¬ 
tions  (A-15)  through  (A-18)  were  used.  Figures  A-2  through 
A-7  show  the  results  of  the  tests. 


136 


Figure  A-2.  cos  2ttx/6 
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Figure  A-5.  (sin  2ux/4)(sin  2uy/24) 


Figure 


Figure  A-l  had  as  its  input  a  24  x  24  array  with  magni¬ 
tudes  of  each  element  varying  only  across  the  rows  via 
cos(2irx/6).  Using  the  above  interpretation  of  a  two-dimen¬ 
sional  Fourier  Transform,  the  fundamental  frequency  assumed 
would  be  one  which  corresponded  to  one  period  across  the  24 
element  array.  Since  the  input  obviously  fits  4  cycles  in 
that  space,  the  only  nonzero  component  of  the  Fourier  Trans¬ 
form  expected  would  be  at  the  fourth  harmonic.  The  transform 
would  also  be  expected  to  be  real  only  since  the  input  is  a 
pure  cosine.  Figure  A-2  is  then  encouraging  in  that  the  only 
nonzero  components  are  where  expected  and  the  imaginary  por¬ 
tions  of  that  answer  is  extremely  small.  The  reason  that 
there  are  nonzero  elements  only  in  the  first  row  is  that  when 
the  transform  was  accomplished  across  the  columns,  there  was 
no  variation  in  the  y  coordinate,  which  is  equivalent  to  taking 
the  transform  of  a  constant  one.  That  transform  results  in  a 
one  in  the  zero  frequency  components  in  each  column  which  is 
the  first  row  and  zeros  everyplace  else.  When  the  multipli¬ 
cation  of  Equation  (A-9)  is  accomplished.  Figure  A-2  results. 

Figure  A-3  had  an  input  of  (cos  2irx/6)  •  (cos  2rry/24) 
which  is  the  same  variation  in  the  x-direction  and  a  varia¬ 
tion  in  the  y-direction  which  corresponds  to  one  period 
exactly.  The  only  nonzero  result  from  the  transform  along 
the  columns,  variation  in  y,  should  be  in  the  fundamental. 


which  it  is 


Figure  A-4  through  A-7  are  just  further  examples  to 
demonstrate  Fourt.  Figure  A-4  had  an  input  of  the  fourth 
harmonic  in  the  x-direction  and  the  fundamental  in  the  y- 
direction.  Figure  A-5  contains  the  sixth  harmonic  in  the 
x-direction.  Figure  A-6  contains  the  twelfth  harmonic  in 
the  x-direction  which  only  appears  as  a  conjugate.  Figure 
A-7  contains  only  the  twelfth  harmonic  in  both  directions 
and  encloses  spatial  frequencies  within  boxes.  A  thorough 
understanding  of  these  results  is  needed  to  understand  taking 
derivatives  and  shifting  in  the  transform  domain. 

In  summary  Figure  A-8  is  provided  to  show  the  output 
of  the  subroutine  Fourt.  The  DC  component  is  the  product  of 
the  zero  frequency  components  of  the  one-dimensional  trans¬ 
forms  in  both  directions.  Elements  Tj1  through  T^24  24^ 
can  be  interpreted  similarly  to  Equations  (A-10)  through 
(A-l  3)  . 


144 


liiHHiniiH 

■■■■■■■■■■■ ■■■ 

SSSSSSSSSBSSa 


■■■■■■I 

iiiiSii 


11  harmonic 
in  y  variation 

x.  o_ 

conjugate  12 
harmonic  in  y 


conjugate  1st 
harmonic  in  y 


24,24) 


>1  c 

C  -C 

o  0 

U  0  -P  X 

C  -H 

■H  -H  <N 

0  4-> 

C  -P  -H  C 

0  id 

0  u  -H 

tr-H 

e  <o  a) 

0)  P 

p  p  .p  u 

P  (0 

tO  -H  (0  H 

4-1  > 

A  T3  01  C 
0  0 

0  X 

«  X  •-»  E 

p 

■P  CP 

0)  c 

C  0  fd 

N  -H 

rH  *h  u  n 

conjugate 


Appendix  B 


Shifting  Property  of  the  Two-Dimensional  Finite 
Discrete  Fourier  Transform 


To  accomplish  interframe  filtering,  the  target's  inten¬ 
sity  profile  must  be  centered  within  each  data  frame.  This 
would  allow  the  successively  centered  frames  to  be  averaged 
to  attenuate  noise.  However,  for  the  reasons  presented  in 
Chapter  I,  the  true  target's  intensity  profile  is  offset  from 
the  center  of  the  f ield-of-view.  This  offset  from  the  center 
of  the  data  array  is  estimated  by  the  Kalman  filter  updated 
states,  x(t±+) .  The  problem  then  becomes  how  to  manipulate 
the  data  array  so  that  the  target's  intensity  profile  is 
centered  within  the  field-of-view.  The  solution  is  to  use 
the  shifting  property  of  the  Fourier  Transform  to  alter  the 
data  to  negate  a  spatial  shift  of  an  amount  calculated  from 
the  filter's  updated  estimate. 

As  stated  in  Chapter  II,  a  shift  of  a  function  in  the 
space  domain  introduces  a  .linear  phase  shift  in  the  frequency 
domain  (8:9).  Equation  (8)  of  Chapter  II  is  now  rewritten 
to  include  the  summations  of  the  Two-Dimensional  Finite  Dis¬ 
crete  Fourier  Transform. 


1  N— 1  N-l  ["  1  [" 

F  g(x-xQ,y-y0)  =  Z^g(x,y)exp  — 1(fxx+fyy)  exp  — ^ (fxx0+fyy0) 


0<y0<N 


0<_Xq  <_N 


( B— 1 ) 


This  equation  can  be  shortened,  as  done  with  Equation 

(8)  of  Chapter  II  by  substituting  G(f  , f  )  for  the  Fourier 

x  y 

Transform  of  the  unshifted  function  g(x,y). 


g(x-x0,y-y0) 


=  G(fx/fy)exp 


-i27 r 
N 


(Exx0+£yy0> 


Q<Vn  <N  0<_xn_<N 


(B— 2) 


The  reason  for  the  qualifier  of  shift  magnitudes  under 
Equations  (B-l)  and  (B-2)  is  that  if  a  shift  were  allowed  to 
be  greater  than  the  period  assumed  by  the  Fourier  Transform, 
see  Appendix  A,  it  could  not  be  distinguished  from  a  shorter 
shift.  This  is  because  of  aliasing.  Therefore,  Xq<N  and 
yQ£N  are  required  for  the  interpretation  of  the  shift  as  a 
rotation  in  two  dimensions  to  be  a  valid  unique  representa¬ 
tion. 


The  implication  of  Equations  (B-l)  and  (B-2)  using 
different  factors  for  different  spatial  frequencies  is  that 
the  results  of  Appendix  A,  Two-Dimensional  Fourier  Trans¬ 
forms,  must  be  used  to  locate  each  of  the  spatial  frequencies 
within  the  transform  domain.  To  implement  the  shift,  each 
point  in  the  transform  domain  must  be  multiplied  by  the  con¬ 
jugate  of  the  phase  shift  induced  by  the  spatial  translation. 


g  (x,  y)  =  F 


'1  G (f  , 
x 


fjexp 


"x  0  yJ  0 


exp 


+  i  2tt 


N 


“  F_1  F[9<*-Vy-y0>]  exp  ±SJi(fxWo) 


(£xx0+fyy0>l 

(B-3) 
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The  centered  intensity  function  g(x,y)  can  be  derived 
from  the  inverse  Fourier  Transform  of  the  Fourier  Transform 
of  the  shifted  function  F[g (x-xQ, y-yQ) ]  via  Equation  (B-3) . 
The  shifted  intensity  function  g(x-XQ,y-yg)  is  the  input 
frame  of  data  in  Figure  1  of  Chapter  I.  The  Two-Dimensional 
Finite  Discrete  Fourier  Transform  of  this  data  is  then  taken, 
which  is  F [g (x-xQ, y-yQ) ]  .  In  Appendix  A,  Figure  A-8,  the 
output  of  a  two-dimensional  transform  is  shown.  To  implement 
a  shift  each  point  of  the  two-dimensional  transformed  data 
array  must  be  multiplied  by  the  conjugate  of  the  translation 
induced  linear  phase  shift  as  shown  in  (B-3) . 

Subroutine  Shift  shifts  a  data  array  by  an  amount  spe¬ 
cified  by  its  parameters.  In  Figure  A-8,  for  example,  point 
T(1  2)  tJie  P0-*-11*"  t*ie  transformed  array  in  the  first 
row  second  column.  Figure  A-8  shows  it  is  the  first  harmonic 
in  the  x-variation  and  the  zero  frequency  in  the  y-variation 

(i.e.,  f  =l,f  =0).  This  point  would  therefore  be  multiplied 
x  y 

by  exp  "xo^  to  imPlement  a  shift  of  Xq,  yQ.  Taking 
another  point,  for  example,  T^4  24)  transformed  data 

array  component  in  the  twenty- fourth  row  twenty-fourth  column 
Figure  A-8  shows  this  point  to  be  the  product  of  the  conju¬ 
gate  of  the  first  harmonics  in  both  directions.  To  implement 
a  shift  this  point  would  be  multiplied  by  exp[-^ — (xQ+y0) ] . 

Figure  B-la  and  B-lb  show  the  results  of  subroutine 
Shift.  The  inputs  to  Shift  are  the  data  array,  the  dimension 


8 


(a)  Original  data 


Figure  B 


1  1 
1  1 


shift  =  2.,  y  shift  =  2. 


1.  Results  of  Subroutine  Shift 


i> 


of  that  data  array,  and  the  amount  in  pixels  of  shift  desired. 
This  routine  was  tested  with  data  representing  a  single 
Gaussian  and  a  multi-hotspot  three  Gaussian  profile. 

In  summary.  Equation  (B— 3 )  shows  how  the  centered  in¬ 
tensity  function  can  be  derived  from  the  Fourier  Transform 
of  the  spatially  translated  intensity  function  F[g (x-XQ,y-yQ) ]. 
The  Two-Dimensional  Finite  Discrete  Fourier  Transform  assumes 
that  this  finite-area  array  represents  one  period  of  a  two- 
dimensional  periodic  sequence.  For  this  reason,  if  a  shift 
is  greater  than  the  period  it  can't  be  distinguished  from  a 
shorter  shift. 
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Appendix  C 

Implementation  of  the  Exponential  Smoothing  Algorithm 

As  explained  in  Chapter  II,  the  target's  intensity 
pattern  is  corrupted  by  noise.  Interframe  smoothing  of 
centered  intensity  patterns  is  accomplished  with  an  exponen¬ 
tial  smoothing  algorithm  defined  by 

y (t)  =  ay(t)  +  (a^l)  y(t-l)  (C-l) 

where 

£(t)  =  current  averaged  data  frame 
y(t)  =  current  data  frame 
£(t-l)  =  previous  result  of  averaging  data 
a  =  smoothing  constant  ’ 0<a<l 

For  smoothing  in  the  frequency  domain  Equation  (C-l)  becomes 
£(t  )  =  aL  (t±)  +  (1-a)  L(ti_1)  (C-2) 

where 

/v 

L  (t . )  =  current  averaged  data  frame  in  frequency  do- 
1  main 

L (t . )  =  transform  of  current  noise  corrupted  data 
frame 

A 

kvt-  t)  =  previous  result  of  averaging  in  frequency 
1_  domain 

a  =  smoothing  constant  0<a<l 

If,  for  example,  the  steady  state  value  of  a  is  to  be 
0.1,  this  algorithm  would  be  implemented  as  follows.  The 
smoothing  constant,  a,  is  varied  for  the  first  ten  frames 
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(a=l/k  k=l, 2,3, . . .10)  until  the  steady  state  value  of  a  =  .1 
is  reached.  Once  steady  state  is  reached  Equation  (C-2) 
becomes : 


i(tL)  =  .1  L(ti)  +  .9  L(ti_1)  (C-3) 

when  the  smoothing  constant,  at,  is  set  to  1/k  for  the  first 
ten  frames,  k  is  equal  to  from  one  to  ten,  the  equations 
below  describe  the  smoothing  algorithm. 


k  =  1 

kO^) 

=  k(tx) 

k  =  2 

£(t2) 

=  1 2  ^ 

+  2^(tl} 

2—^2^ 

+  ^kft^) 

k  =  3 

£(t3) 

’  ^(t3> 

+ 

-  3i(t3> 

=  jL(t3)  +  +  Lft^] 

k  =  4  £(t4)  =  ^L(t4)  +  |-£(t3) 


=  Jb(t4)  +  ~tj(L(t3)  +  L(t2)  +L(t1)] 

=  j[L(t4)  +  k(t3)  +  k(t2)  +  (C-4) 

This  time  averaging  of  the  first  k  frames  continues  until 
the  steady  state  value  of  a  =  .1  and  Equation  (C-3)  is 
reached.  For  the  tenth  frame  the  smoothing  algorithm 
becomes : 

k  =  10  L(t10)  =  XoL(t10>  +  lO^^V  +  •  •  •  +  kC^)] 

=  77T[L(t.n)  +  .  .  .  +  L(t  )] 


k  =  11  L(ti;L)  =  ^(t^  +  (.9)(.l)*  E^L(t,) 

k  =  12  k(t12)  =  10-(t12)  +  k(tn)  +  (.9)  (.9)  (.1)  * 

10 

E  L(t  )  (C-5) 

i=l  1 

For  k>ll  the  exponential  smoothing  equation  can  be  generalized 
to  the  following  equation 

a  2.  k-11  ^ 

k>12  L(tk)  =  ~L(tk)  +  ^E  (.9)*  (.1)  Lk_£ 

+  ( . 9) k~10  (.1)  [L(t-10)  +  .  .  .  rt(t1) ]  (C-6) 


When  k=ll  the  second  term  of  Equation  (C-6) ,  the  sum¬ 
mation  from  £=1  to  k-11,  be,  s  zero  since  the  upper  limit 
on  the  summation  is  less  than  the  lower  limit. 

If  it  is  desired  to  do  this  algorithm  in  the  spatial 
domain  the  same  summations  hold. 


i=ll 


i>ll 


i<10 


h(x,y,t.)  =  .1  h(x,y,t,)  +  E  (.9)  (.1) k) 

1  1  k=l  1 

+  (.9) 1-10  ( .1)  [h(x,y,t1Q)  +  .  .  .  +h(x,y,t1)] 
1  i  i 

h  (x,y,  t . )  =FA[L(t.)3  =  7  E  h  (x,  y,  t.  )  (C-7) 

1  1  1  k=l  K 


In  summary,  an  exponential  smoothing  algorithm  was 
used  to  minimize  the  effects  of  the  corrupting  noise.  Equa¬ 
tions  (C— 4 ) ,  (C-5),  and  (C-6)  express  the  implementation  of 
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this  algorithm  in  the  frequency  domain.  This  same  algorithm 
is  valid  in  the  spatial  domain  as  shown  in  Equation  (C— 7 ) . 


* 


Implementation  of  the  Spatial  Derivatives 
Using  Fourier  Transforms 


The  extended  Kalman  filter  algorithm  requires  the  spa¬ 
tial  derivatives  of  the  two-dimensional  intensity  function 
with  respect  to  each  of  the  states  (see  Chapter  IV) .  The 
derivative  property  of  the  Fourier  Transform  can  be  used  to 
obtain  these  derivatives. 


F 


F 


9h (x, v) 


3x 

3h(x,.y.), 

9X 


=  j2Trfx-  F[h  (x,y)  ] 
=  j27Tfy-  F [h  (x, y)  ] 


(D-l) 


Equation  (D-l)  shows  that  differentiation  in  the  x  direction 
of  the  space  domain  is  equivalent  to  multiplication  by  j2v(fx) 
in  the  frequency  domain  and  similarly  for  the  y  direction. 

The  derivative  of  the  intensity  function  can  be  expressed  in 
terms  of  the  transform  of  that  function. 

To  implement  Equation  (D-l),  the  results  of  Appendix  A 
must  be  used  to  locate  each  of  the  spatial  frequencies  in 
the  transform  domain.  To  obtain  the  spatial  derivative  of 
the  intensity  function,  each  data  point  in  the  transform 
domain  must  be  multiplied  by  j 2n* (corresponding  spatial  fre¬ 
quency  in  the  direction  of  the  desired  derivative) .  In 
Figure  A-8,  for  example  point  T^  ^  is  the  point  of  the 
transformed  array  in  the  first  row  and  second  column.  Figure 
A-8  shows  this  point  to  be  the  product  of  the  first  harmonic 
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in  the  x-variation  and  the  zero  frequency  in  the  y-variation. 
To  obtain  the  spatial  derivative  in  the  x-spatial-direction, 
this  point  would  be  multiplied  by  (j2ir/N).  Taking  another 
point,  say  T^4  24)'  is  the  transformed  data  array 

component  in  the  twenty-fourth  row  and  twenty- fourth  column, 
the  conjugate  of  spatial  frequencies  can  be  demonstrated. 

T ( 2 4  24 j  is  the  product  of  the  conjugate  of  the  first  har¬ 
monics  in  both  x  and  y  directions.  To  obtain  the  spatial 
derivative  in  the  x-direction  this  point  would  be  multiplied 
by  -j2ir/N. 

To  show  how  the  Derivative  algorithm  was  tested  Figures 
D— 1  through  D-3  are  included.  Figure  D-l  shows  the  original 
data  to  be  sinusoids  in  both  spatial  directions,  of  one 
period  in  x  and  two  periods  in  y,  sin  (2ttx/24)  *sin  (2iry/12)  . 
Figure  D-2  shows  the  result  of  using  the  Derivative  Property 
of  the  Fourier  Transform  tc  obtain  the  spatial  derivative  in 
the  x-direction.  The  result  is  a  cosinusoid  of  one  period. 
Similarly,  Figure  D-3  shows  the  spatial  derivative  in  the 
y-direction.  The  result  is  a  cosinusoid  of  two  periods. 

These  figures  use  a  grey-scale  routine  which  attempts  to 
show  increases  in  intensity  by  more  over-prints.  (see  Sub¬ 
routine  Display) .  This  routine  only  reflects  th  ^  absolute 
values  of  intensities. 

In  summary,  Equation  D-l  shows  how  the  spatial  deriva¬ 
tive  can  be  expressed  in  terms  of  the  transform  of  that  in- 
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Figure  D- 
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Figure  D-2.  Spatial  Derivative  in  x 
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Figure  D-3.  Spatial  Derivative  in  y 
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tensity  function.  Within  the  computer  simulation,  subroutine 
Deriv  implements  this  algorithm.  This  appendix  derives  the 
spatial  derivative  of  an  intensity  function.  As  shown  in 
Chapter  IV  the  algorithm  will  require  the  negative  of  this 
result. 
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Appendix  E 

Generation  of  White  Gaussian  Noise  Process 

Throughout  this  research,  it  was  desired  to  generate 
samples  of  a  discrete-time  white  Gaussian  noise  vector  pro¬ 
cess  with  a  mean  of  zero  and  a  given  variance.  This  function 
must  be  simulated  through  the  use  of  pseudorandom  codes  that 
generate  numbers  as  though  they  were  generated  as  samples  of 
a  scalar  random  variable  with  uniform  probability  density 
between  0  and  1.  (see  Figure  E-l) 


Figure  E-l.  Output  Probability  Density  of  Pseudorandom  Code 
The  mean  and  variance  of  this  distribution  is  given  by 

oo  2. 

y,  =  /  A  f  (A)dA  =  /1-AdA  =  ~  (E-l) 

A  _co  x  0  Z 

oo  2 

0*  /  (X"2)2fx(X)dX  =  nX2-A+^JdA=^  (E-2) 

—  oo  0 

If  twelve  independent  calls  are  made  to  such  a  random  number 
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generator  and  the  outputs  of  these  calls,  specific  realiza¬ 


tions  are  summed,  as 

12 

A,  =  £  S,  (E-3) 

J  i=l  1 

the  result  is  a  realization  of  a  random  variable  with  a  mean 
of  6  and  a  variance  of  1  that  is  essentially  Gaussian  (by 
the  Central  Limit  Theorem,  Law  of  Large  Numbers,  etc.).  Six 
is  then  subtracted  from  each  sum  of  realizations,  ft  ,  to  pro¬ 
duce  a  discrete  realization  of  E  .  which  will  have  a  Gaussian 

J 

distribution  of  zero  mean  and  unity  variance 

Cj  =  Qj  -  6  (E— 4) 

Now  the  Cj's  are  arrayed  into  a  vector  of  appropriate  dimen¬ 
sion.  For  example  to  create  a  64  x  1  vector  w^,  64  realiza¬ 
tions  of  E, y  each  being  created  from  independent  calls  to 
the  pseudorandom  code,  would  be  arrayed: 


(E— 5) 


The  vector  process  w^  is  composed  of  independent  scalar 
white  Gaussian  noises  of  mean  zero  and  unit  variance. 


To  create  a  discrete-time  white  Gaussian  noise  vector 

process  w,(t.),  described  by  mean  of  zero  and  covariance  of 
“Td  i 

q  (t . ) *  in  general  non-diagonal,  the  following  equation  can 
— d  r 

be  used. 


(E-6) 


where  ^“denotes  Cholesky  lower  triangular  square  root.  (13:408) 
Reference  13  shows  that  Equation  (E-6)  properly  models  the 
desired  characteristics. 

In  summary,  this  appendix  shows  how  samples  of  discrete¬ 
time  white  Gaussian  noise  vectors  were  generated  from  pseudo¬ 
random  computer  codes. 
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Appendix  F 

Discrete  Representation  of  the  Derivative  of  the 
Intensity  Profile  with  Respect  to  the 
Kalman  Filter  States 

To  present  a  clearer  presentation  of  the  discrete  repre¬ 
sentation  of  the  intensity  profile  and  its  derivative  with 
respect  to  the  Kalman  filter  states,  the  continuous  one¬ 
dimensional  Gaussian  profile  of  Figure  F-la  is  used.  If  the 
discrete  representation  of  this  intensity  function  is  approxi¬ 
mated  as  a  single  value  at  the  center  of  a  pixel  the  Figure 
F-lb  is  appropriate. 

At  time  t.,  assume  that  the  intensity  profile  is  cen¬ 
tered  in  the  f ield-of-view  of  eight  horizontal  pixels,  as  in 
Figure  F-lb.  Although  h  is  actually  the  averaged  value  over 
a  pixel  it  is  approximated  as  the  sampled  value  at  the  center. 
At  time  the  target  has  moved  in  the  direction  of  increas¬ 

ing  state  by  an  amount  Xq,  see  Figure  F-2.  Taking  the  first 
pixel  as  an  example,  h^  decreased  in  magnitude  as  the  state 
increased.  The  magnitudes  of  the  intensity  measurements  for 
all  pixels  located  at  state  values  to  the  left  of  the  profile's 
maximum  decreased  in  intensity  while  those  to  the  right  in¬ 
creased  in  magnitude.  This  is  the  spatial  derivative  rotated 
180  degrees  about  the  vertical  axis.  Figure  F-3  gives  the 
spatial  derivative  of  a  one-dimensional  intensity  profile  as 
in  Figure  F-l,  while  Figure  F-4  gives  the  derivative  with 
respect  to  the  Kalman  filter  states. 
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In  summary,  the  derivative  of  the  intensity  profile 
with  respect  to  the  Kalman  filter  states  is  the  negative  of 
the  spatial  derivative  of  that  intensity  profile  with  respect 
to  coordinate  directions. 


x (pixels) 


One-Dimensional  Gaussian 


Appendix  G 
Monte  Carlo  Study 

For  the  Monte  Carlo  analysis,  many  samples  of  the  error 
process  are  generated  within  the  computer  simulation,  and  the 
error  statistics  are  computed  from  those  samples.  In  this 
research,  twenty  FLIR  data  frames  represent  one  tracking 
history  in  time.  For  each  quantity  of  interest,  errors  are 
computed  before  and  after  processing  a  FLIR  data  frame. 

These  critical  values,  which  the  Kalman  filter  estimates, 
are  compared  to  the  truth  model  output  to  compute  errors,  as 
depicted  in  Figure  G-l  (13:327). 


Figure  G-l.  Generation  of  Error  Samples 
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l 


The  truth  model  enerates  the  measurement  vector  z.  (t , ) 

— t  i 

which  is  a  realization  of  the  measurement  process  (• , • )  at 
time  t.,  i.e.  z.  (t  . ,  w.)  .  The  feedback  loop  only  controls 
the  center  of  the  field-of-view  for  which  the  intensity  mea¬ 
surements  are  generated.  Twenty  such  twenty-frame  time  his¬ 
tories  are  now  processed  and  the  statistics  of  the  errors, 
at  each  time  frame  are  computed  by  averaging  across  the  twenty 
time  histories.  For  example  say  the  difference  between  the 
truth  model  value  for  x  dynamics,  *t(j,  and  the  filter  esti¬ 
mated  value  for  x  dynamics,  x^,  in  the  first  time  history, 

2 

at  frame  one  was  e. .  .  By  keeping  stored  e . .  and  e . . 

•Lixd  1  <1  1  Kxd 


for  each  frame  i  in  a  given  simulation  run,  i=l,  .  .  .  ,  20, 
and  for  the  twenty  time  histories  k=l,  2,  ,  20,  the 

mean  error  and  the  variance  of  the  errors  for  any  frame  can 
be  computed. 

Figure  G-l  assumes  that  the  true  values  of  the  critical 
quantities  yt  are  related  to  the  truth  model  states,  xt»  by  a 
linear  transformation  represented  by  the  matrix  C^.  The 
matrix  C.  is  similarly  explained.  For  every  quantity  of  cri¬ 
tical  interest,  the  appropriate  entries  in  the  1 C*  matrices 
designate  what  combination  of  states  constitute  that  value 
so  the  error  sample  can  be  generated. 

In  summary,  the  object  of  this  Monte  Carlo  study  was 
to  characterize  the  error  process  of  the  algorithm  of  Figure 
1  statistically. 
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Computer  Software 


This  appendix  contains  the  Fortran  source  code  for 
some  of  the  computer  programs  written  in  this  study.  The 
first  program  contained  in  this  appendix  was  written  for  use 
on  the  CDC  Fortran  IV  compiler.  This  program  is  the  completed 
implementation  of  the  algorithm  of  Figure  1.  The  next  program 
was  written  for  use  on  the  MODCOMP  Classic  minicomputer  to 
implement  the  algorithm  of  Chapter  V.  The  last  listing  is 
the  program  which  was  written  to  generate  the  plots  of  Chapter 
VI.  This  software  is  well  commented  to  allow  for  ease  of 
understanding  and  modification. 
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Plots  of  Tracking  Errors 


The  sequence  of  these  plots  is  explained  in  Chapter  VI. 
This  appendix  contains  plots  for  four  settings  of  parameters 
which  control  intensity  function  derivation  and  Kalman  filter 
characteristics.  The  legend  on  each  plot  provides  information 
on  the  target's  spread  parameter,  the  number  of  zeros  which 
are  padded,  the  smoothing  constant,  the  number  of  frequencies 
to  zero,  the  background  noise  strength,  the  standard  deviation 
of  the  discrete  time  noise  driving  the  target  dynamic  states 
of  the  filter  respectively. 
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Extended  Kalman  Filters 


This  appendix  provides  an  overview  of  the  generic  ex¬ 
tended  Kalman  filter  algorithm  (17)  .  The  generic  terms  of 
the  algorithm  are  then  presented  in  terms  of  the  tracking 
example. 

The  system  which  contains  the  states  which  are  to  be 
estimated  is  assumed  to  satisfy  the  nonlinear  stochastic 
differential  equation 

x  ( t )  =  f  [x  (t)  ,u  (t)  ,  t]  +  G(t)  w(t)  (J-l) 

where  the  initial  condition  or.  the  state,  x(t  ),  is  modelled 

o 

A 

as  a  Gaussian  random  vector  with  mean  xq  and  covariance  . 
The  dynamic  driving  noise  w(t)  is  a  zero-mean  white  Gaussian 
noise  process  of  strength  t) . 

E{w(t)  wT(t-t-x)}  =  Q.(t)  6(t)  (J-2) 

The  dynamic  driving  noise,  n(t) ,  is  assumed  to  enter  in  a 
linear  additive  fashion  in  Equation  (J-l) . 

For  this  application,  the  generic  Equation  (J-l)  be¬ 
comes  Equation  (4-4)  where  the  generic  f.  is  replaced  by  a 
constant  F.  The  four  filter  states,  target  position  in  x-y 
coordinates  resulting  from  dynamics  or  atmospherics  (x^  y^ 
x  y  ]  ,  are  multiplied  by  the  constant  filter  plant  matrix 

3  3 

V 

The  generic  discrete-time  measurement  equation  is 
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modelled  as  a  noni inear  function  of  the  states  and  time  plus 
linearly  additive  noise  corruption 


Z. (t±)  =  ll[2S.(ti)  ,  t±]  +  v(ti)  ( J -3 ) 

The  measurement  corruption  noise  v(t^)  is  a  white  zero-mean 
Gaussian  noise  of  strength  R(t^)  and  independent  of  the 
dynamic  driving  noise  w(t^)  and  the  initial  conditions  on 
the  states. 


E{v(ti)  vT(ti)} 


R(ti) 

0 


ti±tj 


(J-4 ) 


In  this  application  the  measurement  vector,  jz(t^)  is 
composed  of  the  64  intensity  measurements  from  the  FLIR. 

The  h[*,*]  function  is  the  expected  intensity  measurements 
as  a  function  of  the  estimated  states.  The  derivation  of 
this  intensity  function  is  the  subject  of  Chapter  II. 

The  propagation  and  measurement  update  equations  for 
the  extended  Kalman  filter  are  presented  in  Chapter  IV. 
These  equations  are  coupled  to  the  state  estimate  relations 

A 

via  the  dependence  of  h[*,*]  on  the  state  estimate  x(t2) . 
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